Skip to content

GitLab

Commit 873552ba authored by Peter Eastman's avatar Peter Eastman
Browse files

Deleted free energy plugin

parent c775bd19
Hide whitespace changes
Inline Side-by-side
Showing with 0 additions and 5747 deletions
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateObcGbsaSoftcoreForces2.cu deleted 100755 → 0
View file @ c775bd19
/* -------------------------------------------------------------------------- *
* 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/>. *
* -------------------------------------------------------------------------- */
#include "openmm/OpenMMException.h"
#include "gputypes.h"
#include "freeEnergyGpuTypes.h"
#include <stdio.h>
#include <cuda.h>
#include <vector_functions.h>
#include <cstdlib>
#include <string>
#include <iostream>
#include <fstream>
using namespace std;
struct Atom {
float x;
float y;
float z;
float r;
float sr;
float npScale;
float fx;
float fy;
float fz;
float fb;
};
static __constant__ cudaGmxSimulation cSim;
static __constant__ cudaFreeEnergyGmxSimulation feSimDev;
extern "C"
void SetCalculateObcGbsaSoftcoreForces2Sim( freeEnergyGpuContext freeEnergyGpu )
{
cudaError_t status;
status = cudaMemcpyToSymbol(cSim, &freeEnergyGpu->gpuContext->sim, sizeof(cudaGmxSimulation));
RTERROR(status, "cudaMemcpyToSymbol: SetCalculateObcGbsaSoftcoreForces2Sim copy to cSim failed");
status = cudaMemcpyToSymbol( feSimDev, &freeEnergyGpu->freeEnergySim, sizeof(cudaFreeEnergyGmxSimulation));
RTERROR(status, "cudaMemcpyToSymbol: SetCalculateObcGbsaSoftcoreForces2Sim copy to feSimDev failed");
}
// Include versions of the kernels for N^2 calculations.
#define METHOD_NAME(a, b) a##N2##b
#include "kCalculateObcGbsaSoftcoreForces2.h"
#define USE_OUTPUT_BUFFER_PER_WARP
#undef METHOD_NAME
#define METHOD_NAME(a, b) a##N2ByWarp##b
#include "kCalculateObcGbsaSoftcoreForces2.h"
// Include versions of the kernels with cutoffs.
#undef METHOD_NAME
#undef USE_OUTPUT_BUFFER_PER_WARP
#define USE_CUTOFF
#define METHOD_NAME(a, b) a##Cutoff##b
#include "kCalculateObcGbsaSoftcoreForces2.h"
#define USE_OUTPUT_BUFFER_PER_WARP
#undef METHOD_NAME
#define METHOD_NAME(a, b) a##CutoffByWarp##b
#include "kCalculateObcGbsaSoftcoreForces2.h"
// Include versions of the kernels with periodic boundary conditions.
#undef METHOD_NAME
#undef USE_OUTPUT_BUFFER_PER_WARP
#define USE_PERIODIC
#define METHOD_NAME(a, b) a##Periodic##b
#include "kCalculateObcGbsaSoftcoreForces2.h"
#define USE_OUTPUT_BUFFER_PER_WARP
#undef METHOD_NAME
#define METHOD_NAME(a, b) a##PeriodicByWarp##b
#include "kCalculateObcGbsaSoftcoreForces2.h"
void kCalculateObcGbsaSoftcoreForces2( freeEnergyGpuContext freeEnergyGpu )
{
unsigned int threadsPerBlock;
static unsigned int threadsPerBlockPerMethod[3] = { 0, 0, 0 };
static unsigned int natoms[3] = { 0, 0, 0 };
gpuContext gpu = freeEnergyGpu->gpuContext;
unsigned int methodIndex = static_cast<unsigned int>(freeEnergyGpu->freeEnergySim.nonbondedMethod);
if( methodIndex > 2 ){
throw OpenMM::OpenMMException( "kCalculateObcGbsaSoftcoreForces2 method index invalid." );
}
if( natoms[methodIndex] != gpu->natoms ){
unsigned int extra = methodIndex == 0 ? 0 : sizeof(float3);
threadsPerBlockPerMethod[methodIndex] = std::min(getThreadsPerBlockFEP( freeEnergyGpu, (sizeof(Atom) + extra), gpu->sharedMemoryPerBlock ), gpu->sim.nonbond_threads_per_block );
natoms[methodIndex] = gpu->natoms;
}
threadsPerBlock = threadsPerBlockPerMethod[methodIndex];
switch (freeEnergyGpu->freeEnergySim.nonbondedMethod)
{
case FREE_ENERGY_NO_CUTOFF:
if (gpu->bOutputBufferPerWarp)
kCalculateObcGbsaSoftcoreN2ByWarpForces2_kernel<<<gpu->sim.bornForce2_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit);
else
kCalculateObcGbsaSoftcoreN2Forces2_kernel<<<gpu->sim.bornForce2_blocks, threadsPerBlock,
sizeof(Atom)*threadsPerBlock>>>(gpu->sim.pWorkUnit);
break;
case FREE_ENERGY_CUTOFF:
if (gpu->bOutputBufferPerWarp)
kCalculateObcGbsaSoftcoreCutoffByWarpForces2_kernel<<<gpu->sim.bornForce2_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float3))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit);
else
kCalculateObcGbsaSoftcoreCutoffForces2_kernel<<<gpu->sim.bornForce2_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float3))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit);
break;
case FREE_ENERGY_PERIODIC:
if (gpu->bOutputBufferPerWarp)
kCalculateObcGbsaSoftcorePeriodicByWarpForces2_kernel<<<gpu->sim.bornForce2_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float3))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit);
else
kCalculateObcGbsaSoftcorePeriodicForces2_kernel<<<gpu->sim.bornForce2_blocks, threadsPerBlock,
(sizeof(Atom)+sizeof(float3))*threadsPerBlock>>>(gpu->sim.pInteractingWorkUnit);
break;
}
LAUNCHERROR("kCalculateObcGbsaSoftcoreForces2");
}
plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateObcGbsaSoftcoreForces2.h deleted 100644 → 0
View file @ c775bd19
/* -------------------------------------------------------------------------- *
* 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 the second stage of GBSA. It is included
* several times in kCalculateObcGbsaSoftcoreForces2.cu with different #defines to generate
* different versions of the kernels.
*/
__global__
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_BORNFORCE2_THREADS_PER_BLOCK, 1)
#elif (__CUDA_ARCH__ >= 120)
__launch_bounds__(GT2XX_BORNFORCE2_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_BORNFORCE2_THREADS_PER_BLOCK, 1)
#endif
void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned int* workUnit)
{
extern __shared__ Atom sA[];
unsigned int totalWarps = gridDim.x*blockDim.x/GRID;
unsigned int warp = (blockIdx.x*blockDim.x+threadIdx.x)/GRID;
unsigned int numWorkUnits = cSim.pInteractionCount[0];
unsigned int pos = warp*numWorkUnits/totalWarps;
unsigned int end = (warp+1)*numWorkUnits/totalWarps;
#ifdef USE_CUTOFF
float3* tempBuffer = (float3*) &sA[blockDim.x];
#endif
unsigned int lasty = -0xFFFFFFFF;
while (pos < end)
{
// Extract cell coordinates from appropriate work unit
unsigned int x = workUnit[pos];
unsigned int y = ((x >> 2) & 0x7fff) << GRIDBITS;
x = (x >> 17) << GRIDBITS;
unsigned int tgx = threadIdx.x & (GRID - 1);
unsigned int i = x + tgx;
float4 apos = cSim.pPosq[i];
float2 a = cSim.pObcData[i];
float fb = cSim.pBornForce[i];
float nonPolarScaleDataI = feSimDev.pNonPolarScalingFactors[i];
unsigned int tbx = threadIdx.x - tgx;
unsigned int tj = tgx;
Atom* psA = &sA[tbx];
sA[threadIdx.x].fx = 0.0f;
sA[threadIdx.x].fy = 0.0f;
sA[threadIdx.x].fz = 0.0f;
float3 af;
af.x = 0.0f;
af.y = 0.0f;
af.z = 0.0f;
if (x == y) // Handle diagonals uniquely at 50% efficiency
{
// Read fixed atom data into registers and GRF
sA[threadIdx.x].x = apos.x;
sA[threadIdx.x].y = apos.y;
sA[threadIdx.x].z = apos.z;
sA[threadIdx.x].r = a.x;
sA[threadIdx.x].sr = a.y;
sA[threadIdx.x].fb = fb;
sA[threadIdx.x].npScale = nonPolarScaleDataI;
for (unsigned int j = (tgx+1)&(GRID-1); j != tgx; j = (j+1)&(GRID-1))
{
float dx = psA[j].x - apos.x;
float dy = psA[j].y - apos.y;
float dz = psA[j].z - apos.z;
#ifdef USE_PERIODIC
dx -= floorf(dx/cSim.periodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
dy -= floorf(dy/cSim.periodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
dz -= floorf(dz/cSim.periodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
#endif
float r2 = dx * dx + dy * dy + dz * dz;
float r = sqrtf(r2);
// Atom I Born forces and sum
float rScaledRadiusJ = r + psA[j].sr;
float l_ij = 1.0f / max(a.x, fabs(r - psA[j].sr));
float u_ij = 1.0f / rScaledRadiusJ;
float rInverse = 1.0f / r;
float l_ij2 = l_ij * l_ij;
float u_ij2 = u_ij * u_ij;
float r2Inverse = rInverse * rInverse;
float t1 = logf(u_ij / l_ij);
float t2 = (l_ij2 - u_ij2);
float t3 = t2 * rInverse;
t1 *= rInverse;
// Born Forces term
float term = 0.125f * (1.000f + psA[j].sr * psA[j].sr * r2Inverse) * t3 +
0.250f * t1 * r2Inverse;
term *= psA[j].npScale;
float dE = fb * term;
#if defined USE_CUTOFF
if (a.x >= rScaledRadiusJ || i >= cSim.atoms || x+j >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)
#else
if (a.x >= rScaledRadiusJ || i >= cSim.atoms || x+j >= cSim.atoms )
#endif
{
dE = 0.0f;
}
float d = dx * dE;
af.x -= d;
psA[j].fx += d;
d = dy * dE;
af.y -= d;
psA[j].fy += d;
d = dz * dE;
af.z -= d;
psA[j].fz += d;
}
// Write results
float4 of;
#ifdef USE_OUTPUT_BUFFER_PER_WARP
unsigned int offset = x + tgx + warp*cSim.stride;
#else
unsigned int offset = x + tgx + (x >> GRIDBITS) * cSim.stride;
#endif
of = cSim.pForce4[offset];
of.x += af.x + sA[threadIdx.x].fx;
of.y += af.y + sA[threadIdx.x].fy;
of.z += af.z + sA[threadIdx.x].fz;
of.w = 0.0f;
cSim.pForce4[offset] = of;
} else {
// Read fixed atom data into registers and GRF
if (lasty != y)
{
unsigned int j = y + tgx;
float4 temp = cSim.pPosq[j];
float2 temp1 = cSim.pObcData[j];
sA[threadIdx.x].fb = cSim.pBornForce[j];
sA[threadIdx.x].npScale = feSimDev.pNonPolarScalingFactors[j];
sA[threadIdx.x].x = temp.x;
sA[threadIdx.x].y = temp.y;
sA[threadIdx.x].z = temp.z;
sA[threadIdx.x].r = temp1.x;
sA[threadIdx.x].sr = temp1.y;
}
float sr2 = a.y * a.y;
#ifdef USE_CUTOFF
unsigned int flags = cSim.pInteractionFlag[pos];
if (flags == 0)
{
// No interactions in this block.
}
//else if (flags == 0xFFFFFFFF)
else if (flags)
#endif
{
// Compute all interactions within this block.
for (unsigned int j = 0; j < GRID; j++)
{
float dx = psA[tj].x - apos.x;
float dy = psA[tj].y - apos.y;
float dz = psA[tj].z - apos.z;
#ifdef USE_PERIODIC
dx -= floorf(dx/cSim.periodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
dy -= floorf(dy/cSim.periodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
dz -= floorf(dz/cSim.periodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
#endif
float r2 = dx * dx + dy * dy + dz * dz;
float r = sqrtf(r2);
// Interleaved Atom I and J Born Forces and sum components
float r2Inverse = 1.0f / r2;
float rScaledRadiusJ = r + psA[tj].sr;
float rScaledRadiusI = r + a.y;
float rInverse = 1.0f / r;
float l_ijJ = 1.0f / max(a.x, fabs(r - psA[tj].sr));
float l_ijI = 1.0f / max(psA[tj].r, fabs(r - a.y));
float u_ijJ = 1.0f / rScaledRadiusJ;
float u_ijI = 1.0f / rScaledRadiusI;
float l_ij2J = l_ijJ * l_ijJ;
float l_ij2I = l_ijI * l_ijI;
float u_ij2J = u_ijJ * u_ijJ;
float u_ij2I = u_ijI * u_ijI;
float t1J = logf(u_ijJ / l_ijJ);
float t1I = logf(u_ijI / l_ijI);
float t2J = (l_ij2J - u_ij2J);
float t2I = (l_ij2I - u_ij2I);
float t3J = t2J * rInverse;
float t3I = t2I * rInverse;
t1J *= rInverse;
t1I *= rInverse;
// Born Forces term
float term = 0.125f *
(1.000f + psA[tj].sr * psA[tj].sr * r2Inverse) * t3J +
0.250f * t1J * r2Inverse;
term *= psA[tj].npScale;
float dE = fb * term;
#if defined USE_CUTOFF
if (a.x >= rScaledRadiusJ || i >= cSim.atoms || y+tj >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)
#else
if (a.x >= rScaledRadiusJ || i >= cSim.atoms || y+tj >= cSim.atoms)
#endif
{
dE = 0.0f;
}
float d = dx * dE;
af.x -= d;
psA[tj].fx += d;
d = dy * dE;
af.y -= d;
psA[tj].fy += d;
d = dz * dE;
af.z -= d;
psA[tj].fz += d;
// Atom J Born sum term
term = 0.125f *
(1.000f + sr2 * r2Inverse) * t3I +
0.250f * t1I * r2Inverse;
term *= nonPolarScaleDataI;
dE = psA[tj].fb * term;
float rj = psA[tj].r;
#ifdef USE_CUTOFF
if (rj >= rScaledRadiusI || i >= cSim.atoms || y+tj >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)
#else
if (rj >= rScaledRadiusI || i >= cSim.atoms || y+tj >= cSim.atoms )
#endif
{
dE = 0.0f;
}
dx *= dE;
dy *= dE;
dz *= dE;
psA[tj].fx += dx;
psA[tj].fy += dy;
psA[tj].fz += dz;
af.x -= dx;
af.y -= dy;
af.z -= dz;
tj = (tj + 1) & (GRID - 1);
}
}
#ifdef USE_CUTOFF
else
{
// Compute only a subset of the interactions in this block.
for (unsigned int j = 0; j < GRID; j++)
{
if ((flags&(1<<j)) != 0)
{
float dx = psA[j].x - apos.x;
float dy = psA[j].y - apos.y;
float dz = psA[j].z - apos.z;
#ifdef USE_PERIODIC
dx -= floorf(dx/cSim.periodicBoxSizeX+0.5f)*cSim.periodicBoxSizeX;
dy -= floorf(dy/cSim.periodicBoxSizeY+0.5f)*cSim.periodicBoxSizeY;
dz -= floorf(dz/cSim.periodicBoxSizeZ+0.5f)*cSim.periodicBoxSizeZ;
#endif
float r2 = dx * dx + dy * dy + dz * dz;
float r = sqrtf(r2);
// Interleaved Atom I and J Born Forces and sum components
float r2Inverse = 1.0f / r2;
float rScaledRadiusJ = r + psA[j].sr;
float rScaledRadiusI = r + a.y;
float rInverse = 1.0f / r;
float l_ijJ = 1.0f / max(a.x, fabs(r - psA[j].sr));
float l_ijI = 1.0f / max(psA[j].r, fabs(r - a.y));
float u_ijJ = 1.0f / rScaledRadiusJ;
float u_ijI = 1.0f / rScaledRadiusI;
float l_ij2J = l_ijJ * l_ijJ;
float l_ij2I = l_ijI * l_ijI;
float u_ij2J = u_ijJ * u_ijJ;
float u_ij2I = u_ijI * u_ijI;
float t1J = logf(u_ijJ / l_ijJ);
float t1I = logf(u_ijI / l_ijI);
float t2J = (l_ij2J - u_ij2J);
float t2I = (l_ij2I - u_ij2I);
float t3J = t2J * rInverse;
float t3I = t2I * rInverse;
t1J *= rInverse;
t1I *= rInverse;
// Born Forces term
float term = 0.125f *
(1.000f + psA[j].sr * psA[j].sr * r2Inverse) * t3J +
0.250f * t1J * r2Inverse;
term *= psA[j].npScale;
float dE = fb * term;
#if defined USE_PERIODIC
if (a.x >= rScaledRadiusJ || i >= cSim.atoms || y+j >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)
#elif defined USE_CUTOFF
if (a.x >= rScaledRadiusJ || r2 > cSim.nonbondedCutoffSqr)
#else
if (a.x >= rScaledRadiusJ)
#endif
{
dE = 0.0f;
}
float d = dx * dE;
af.x -= d;
tempBuffer[threadIdx.x].x = d;
d = dy * dE;
af.y -= d;
tempBuffer[threadIdx.x].y = d;
d = dz * dE;
af.z -= d;
tempBuffer[threadIdx.x].z = d;
// Atom J Born sum term
term = 0.125f *
(1.000f + sr2 * r2Inverse) * t3I +
0.250f * t1I * r2Inverse;
term *= nonPolarScaleDataI;
dE = psA[j].fb * term;
float rj = psA[j].r;
#ifdef USE_PERIODIC
if (rj >= rScaledRadiusI || i >= cSim.atoms || y+j >= cSim.atoms || r2 > cSim.nonbondedCutoffSqr)
#elif defined USE_CUTOFF
if (rj >= rScaledRadiusI || r2 > cSim.nonbondedCutoffSqr)
#else
if (rj >= rScaledRadiusI)
#endif
{
dE = 0.0f;
}
dx *= dE;
dy *= dE;
dz *= dE;
tempBuffer[threadIdx.x].x += dx;
tempBuffer[threadIdx.x].y += dy;
tempBuffer[threadIdx.x].z += dz;
af.x -= dx;
af.y -= dy;
af.z -= dz;
// Sum the forces on atom j.
if (tgx % 2 == 0)
{
tempBuffer[threadIdx.x].x += tempBuffer[threadIdx.x+1].x;
tempBuffer[threadIdx.x].y += tempBuffer[threadIdx.x+1].y;
tempBuffer[threadIdx.x].z += tempBuffer[threadIdx.x+1].z;
}
if (tgx % 4 == 0)
{
tempBuffer[threadIdx.x].x += tempBuffer[threadIdx.x+2].x;
tempBuffer[threadIdx.x].y += tempBuffer[threadIdx.x+2].y;
tempBuffer[threadIdx.x].z += tempBuffer[threadIdx.x+2].z;
}
if (tgx % 8 == 0)
{
tempBuffer[threadIdx.x].x += tempBuffer[threadIdx.x+4].x;
tempBuffer[threadIdx.x].y += tempBuffer[threadIdx.x+4].y;
tempBuffer[threadIdx.x].z += tempBuffer[threadIdx.x+4].z;
}
if (tgx % 16 == 0)
{
tempBuffer[threadIdx.x].x += tempBuffer[threadIdx.x+8].x;
tempBuffer[threadIdx.x].y += tempBuffer[threadIdx.x+8].y;
tempBuffer[threadIdx.x].z += tempBuffer[threadIdx.x+8].z;
}
if (tgx == 0)
{
psA[j].fx += tempBuffer[threadIdx.x].x + tempBuffer[threadIdx.x+16].x;
psA[j].fy += tempBuffer[threadIdx.x].y + tempBuffer[threadIdx.x+16].y;
psA[j].fz += tempBuffer[threadIdx.x].z + tempBuffer[threadIdx.x+16].z;
}
}
}
}
#endif
// Write results
float4 of;
#ifdef USE_OUTPUT_BUFFER_PER_WARP
unsigned int offset = x + tgx + warp*cSim.stride;
#else
unsigned int offset = x + tgx + (y >> GRIDBITS) * cSim.stride;
#endif
of = cSim.pForce4[offset];
of.x += af.x;
of.y += af.y;
of.z += af.z;
cSim.pForce4[offset] = of;
#ifdef USE_OUTPUT_BUFFER_PER_WARP
offset = y + tgx + warp*cSim.stride;
#else
offset = y + tgx + (x >> GRIDBITS) * cSim.stride;
#endif
of = cSim.pForce4[offset];
of.x += sA[threadIdx.x].fx;
of.y += sA[threadIdx.x].fy;
of.z += sA[threadIdx.x].fz;
cSim.pForce4[offset] = of;
}
lasty = y;
pos++;
}
}
plugins/freeEnergy/platforms/cuda/src/kernels/kSoftcoreLJ.h deleted 100644 → 0
View file @ c775bd19
#ifndef _K_SOFTCORE_LJ__H__
#define _K_SOFTCORE_LJ__H__
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 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 kernel for calculating softcore LJ force prefactor
*/
#ifdef USE_SOFTCORE_LJ
static __device__ float getSoftCoreLJ( float r2, float sig, float eps, float lambdaI, float lambdaJ, float* energy)
{
float lambda = lambdaI < lambdaJ ? lambdaI : lambdaJ;
eps *= lambda;
// (r/sig)
float sig2 = 1.0f/sig;
sig2 *= sig2;
sig2 *= r2;
float sig6 = sig2*sig2*sig2;
float softcoreLJTerm = 0.5f*( 1.0f - lambda) + sig6;
float softcoreLJInv = 1.0f/softcoreLJTerm;
float softcoreLJInv2 = softcoreLJInv*softcoreLJInv;
*energy = eps*(softcoreLJInv2 - softcoreLJInv);
return eps*softcoreLJInv2*( 12.0f*softcoreLJInv - 6.0f )*sig6;
}
static __device__ float getSoftCoreLJMod( float sigInvR, float eps, float lambdaI, float lambdaJ, float* energy)
{
float lambda = lambdaI < lambdaJ ? lambdaI : lambdaJ;
eps *= lambda;
// (r/sig)
float sig2 = sigInvR*sigInvR;
float sig6 = sig2*sig2*sig2;
float softcoreLJTerm = 0.5f*( 1.0f - lambda) + sig6;
float softcoreLJInv = 1.0f/softcoreLJTerm;
float softcoreLJInv2 = softcoreLJInv*softcoreLJInv;
*energy = eps*(softcoreLJInv2 - softcoreLJInv);
return eps*softcoreLJInv2*( 12.0f*softcoreLJInv - 6.0f )*sig6;
}
#endif
#endif
plugins/freeEnergy/platforms/cuda/staticTarget/CMakeLists.txt deleted 100644 → 0
View file @ c775bd19
#
# Include CUDA related files.
#
SET(OPENMM_BUILD_FREE_ENERGY_PATH ${CMAKE_SOURCE_DIR}/plugins/freeEnergy)
SET(CUDA_NVCC_BUILD_FLAGS "-DCUDPP_STATIC_LIB" )
# INCLUDE(${FINDCUDA_DIR}/FindCuda.cmake)
INCLUDE_DIRECTORIES(${CUDA_INCLUDE})
LINK_DIRECTORIES(${CUDA_TARGET_LINK})
FOREACH(subdir ${OPENMM_SOURCE_SUBDIRS})
FILE(GLOB src_files ${OPENMM_BUILD_FREE_ENERGY_PATH}/platforms/cuda/${subdir}/src/*.cu ${OPENMM_BUILD_FREE_ENERGY_PATH}/platforms/cuda/${subdir}/src/*/*.cu)
SET(SOURCE_FILES ${SOURCE_FILES} ${src_files})
CUDA_INCLUDE_DIRECTORIES(BEFORE ${OPENMM_BUILD_FREE_ENERGY_PATH}/platforms/cuda/${subdir}/include)
CUDA_INCLUDE_DIRECTORIES(BEFORE ${OPENMM_BUILD_FREE_ENERGY_PATH}/platforms/cuda/${subdir}/src)
ENDFOREACH(subdir)
#CUDA_INCLUDE_DIRECTORIES(BEFORE ${OPENMM_BUILD_FREE_ENERGY_PATH}/jama/include)
CUDA_INCLUDE_DIRECTORIES(${OPENMM_BUILD_FREE_ENERGY_PATH}/platforms/cuda/../src
${OPENMM_DIR}/platforms/cuda/src
${OPENMM_DIR}/platforms/cuda/include
${OPENMM_DIR}/platforms/cuda/src/kernels
${OPENMM_DIR}/openmmapi/include )
CUDA_ADD_LIBRARY(${STATIC_TARGET} STATIC ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ${API_ABS_INCLUDE_FILES})
# required for getting OPENMM_EXPORT to be set correctly in 'class OPENMM_EXPORT CudaStreamFactory', ...
# see OpenMM/openmmapi/include/internal/windowsExport.h for details
SET(CUDA_STATIC_COMPILE_FLAG "-DOPENMMCUDAFREEENERGY_BUILDING_STATIC_LIBRARY -DOPENMMCUDAFREEENERGY_USE_STATIC_LIBRARIES -DCUDPP_STATIC_LIB")
SET_TARGET_PROPERTIES(${STATIC_TARGET} PROPERTIES COMPILE_FLAGS ${CUDA_STATIC_COMPILE_FLAG})
TARGET_LINK_LIBRARIES(${STATIC_TARGET} optimized ${OPENMM_LIBRARY_NAME}_static)
TARGET_LINK_LIBRARIES(${STATIC_TARGET} debug ${OPENMM_LIBRARY_NAME}_static_d optimized ${OPENMM_LIBRARY_NAME}_static ${CUFFT_TARGET_LINK})
INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${STATIC_TARGET})
plugins/freeEnergy/platforms/cuda/tests/CMakeLists.txt deleted 100644 → 0
View file @ c775bd19
#
# Testing
#
ENABLE_TESTING()
# INCLUDE(${CMAKE_SOURCE_DIR}/platforms/cuda/cuda-cmake/FindCuda.cmake)
INCLUDE_DIRECTORIES(${CUDA_INCLUDE})
INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/include)
INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/src)
INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/src/kernels)
# serialize test cases for GBVI and GBSAOBC softcore runs
# if INCLUDE_SERIALIZATION is TRUE
SET( INCLUDE_SERIALIZATION FALSE )
#SET( INCLUDE_SERIALIZATION TRUE )
SET( SHARED_OPENMM_TARGET OpenMMFreeEnergy)
SET( SHARED_CUDA_TARGET OpenMMCuda)
IF( INCLUDE_SERIALIZATION )
INCLUDE_DIRECTORIES(${OPENMM_DIR}/serialization/include)
SET( SHARED_OPENMM_SERIALIZATION "OpenMMSerialization" )
SET( SHARED_FREE_ENERGY_SERIALIZATION "FreeEnergySerialization" )
ENDIF( INCLUDE_SERIALIZATION )
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(SHARED_CUDA_TARGET ${SHARED_CUDA_TARGET}_d)
SET(SHARED_OPENMM_TARGET ${SHARED_OPENMM_TARGET}_d)
IF( INCLUDE_SERIALIZATION )
SET(SHARED_OPENMM_SERIALIZATION ${SHARED_OPENMM_SERIALIZATION}_d)
SET(SHARED_FREE_ENERGY_SERIALIZATION ${SHARED_FREE_ENERGY_SERIALIZATION}_d)
ENDIF( INCLUDE_SERIALIZATION )
ENDIF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
# Automatically create tests using files named "Test*.cpp"
FILE(GLOB TEST_PROGS "*Test*.cpp")
FOREACH(TEST_PROG ${TEST_PROGS})
GET_FILENAME_COMPONENT(TEST_ROOT ${TEST_PROG} NAME_WE)
# Link with shared library
CUDA_ADD_EXECUTABLE(${TEST_ROOT} ${TEST_PROG})
IF( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_CUDA_TARGET} ${SHARED_OPENMM_SERIALIZATION} ${SHARED_FREE_ENERGY_SERIALIZATION})
ELSE( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${TEST_ROOT} ${SHARED_TARGET} ${SHARED_OPENMM_TARGET} ${SHARED_CUDA_TARGET})
ENDIF( INCLUDE_SERIALIZATION )
SET(DEFINE_STRING "-DUSE_SOFTCORE")
IF( INCLUDE_SERIALIZATION )
SET(DEFINE_STRING "${DEFINE_STRING} -DOPENMM_SERIALIZE")
ENDIF( INCLUDE_SERIALIZATION )
IF( ${TEST_ROOT} STREQUAL "TestCudaGBVISoftcoreForce2" )
# serialize
SET(DEFINE_STRING "-DTEST_PLATFORM=0 ")
IF( INCLUDE_SERIALIZATION )
SET(DEFINE_STRING "${DEFINE_STRING} -DOPENMM_SERIALIZE ")
SET(TARGET_LINK_LIBRARIES_STRING "${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION}")
ELSE( INCLUDE_SERIALIZATION )
SET(TARGET_LINK_LIBRARIES_STRING "${SHARED_TARGET}")
ENDIF( INCLUDE_SERIALIZATION )
# obc
SET(OBC_DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=1")
SET(OBC_TEST "TestCudaGBSAOBCSoftcoreForce2")
CUDA_ADD_EXECUTABLE(${OBC_TEST} ${TEST_PROG})
IF( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${OBC_TEST} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
ELSE( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${OBC_TEST} ${SHARED_TARGET})
ENDIF( INCLUDE_SERIALIZATION )
SET_TARGET_PROPERTIES(${OBC_TEST} PROPERTIES COMPILE_FLAGS ${OBC_DEFINE_STRING} )
ADD_TEST(${OBC_TEST} ${EXECUTABLE_OUTPUT_PATH}/${OBC_TEST})
# nonbond
SET(NONBOND_DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=0")
SET(NONBOND_TEST "TestCudaNonbondSoftcoreForce2")
CUDA_ADD_EXECUTABLE(${NONBOND_TEST} ${TEST_PROG})
IF( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET} ${SHARED_OPENMM_SERIALIZATION} )
ELSE( INCLUDE_SERIALIZATION )
TARGET_LINK_LIBRARIES(${NONBOND_TEST} ${SHARED_TARGET})
ENDIF( INCLUDE_SERIALIZATION )
SET_TARGET_PROPERTIES(${NONBOND_TEST} PROPERTIES COMPILE_FLAGS ${NONBOND_DEFINE_STRING} )
ADD_TEST(${NONBOND_TEST} ${EXECUTABLE_OUTPUT_PATH}/${NONBOND_TEST})
# gbvi
SET(DEFINE_STRING "${DEFINE_STRING} -DIMPLICIT_SOLVENT=2")
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES COMPILE_FLAGS ${DEFINE_STRING} )
ENDIF( ${TEST_ROOT} STREQUAL "TestCudaGBVISoftcoreForce2" )
#MESSAGE( "${TEST_ROOT} ${DEFINE_STRING}" )
SET_TARGET_PROPERTIES(${TEST_ROOT} PROPERTIES COMPILE_FLAGS ${DEFINE_STRING} )
ADD_TEST(${TEST_ROOT} ${EXECUTABLE_OUTPUT_PATH}/${TEST_ROOT})
ENDFOREACH(TEST_PROG ${TEST_PROGS})
plugins/freeEnergy/platforms/cuda/tests/TestCudaGBVISoftcoreForce.cpp deleted 100644 → 0
View file @ c775bd19
/* -------------------------------------------------------------------------- *
* 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-2009 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 reference implementation of GBVIForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/System.h"
#include "openmm/Context.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "OpenMM.h"
#include "openmm/GBVIForce.h"
#include "openmm/GBVISoftcoreForce.h"
#include "openmm/NonbondedSoftcoreForce.h"
#include "OpenMMFreeEnergy.h"
#include "openmm/freeEnergyKernels.h"
#include "sfmt/SFMT.h"
#include "openmm/VerletIntegrator.h"
#include <iostream>
#include <vector>
#include <algorithm>
#include <cstdio>
#include <typeinfo>
#include <iomanip>
typedef std::vector<double> DoubleVector;
typedef DoubleVector::iterator DoubleVectorI;
typedef DoubleVector::const_iterator DoubleVectorCI;
typedef std::vector<DoubleVector> DoubleVectorVector;
typedef std::pair<int, double> IntDoublePair;
typedef std::vector<IntDoublePair> IntDoublePairVector;
typedef IntDoublePairVector::iterator IntDoublePairVectorI;
typedef IntDoublePairVector::const_iterator IntDoublePairVectorCI;
extern "C" void registerFreeEnergyCudaKernelFactories();
using namespace OpenMM;
using namespace std;
void testSingleParticle( FILE* log ) {
System system;
system.addParticle(2.0);
VerletIntegrator integrator(0.01);
GBVISoftcoreForce* forceField = new GBVISoftcoreForce;
double charge = 1.0;
double radius = 0.15;
double gamma = 1.0;
forceField->addParticle(charge, radius, gamma);
system.addForce(forceField);
NonbondedSoftcoreForce* nonbonded = new NonbondedSoftcoreForce();
nonbonded->setNonbondedMethod(NonbondedSoftcoreForce::NoCutoff);
nonbonded->addParticle( charge, 1.0, 0.0);
system.addForce(nonbonded);
Context context(system, integrator, Platform::getPlatformByName( "Cuda") );
vector<Vec3> positions(1);
positions[0] = Vec3(0, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double bornRadius = radius;
double eps0 = EPSILON0;
double tau = (1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric());
double bornEnergy = (-charge*charge/(8*PI_M*eps0))*tau/bornRadius;
double nonpolarEnergy = -gamma*tau*std::pow( radius/bornRadius, 3.0);
double expectedE = (bornEnergy+nonpolarEnergy);
double obtainedE = state.getPotentialEnergy();
double diff = fabs( obtainedE - expectedE );
if( log ){
(void) fprintf( log, "testSingleParticle expected=%14.6e obtained=%14.6e diff=%14.6e breakdown:[%14.6e %14.6e]\n",
expectedE, obtainedE, diff, bornEnergy, nonpolarEnergy );
}
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
}
/**
* Predicate for sorting <int,double> pair
*
* @param d1 first IntDoublePair to compare
* @param d2 second IntDoublePair to compare
*
*/
bool TestIntDoublePair( const IntDoublePair& d1, const IntDoublePair& d2 ){
return d1.second < d2.second;
}
/**---------------------------------------------------------------------------------------
*
* Get relative difference between two forces
*
* @param f1 force1
* @param f2 force2
* @param forceNorm1 output norm of force1
* @param forceNorm2 output norm of force2
* @param relativeDiff output relative difference between force norms
* @param log if set, output forces
*
*
--------------------------------------------------------------------------------------- */
static void getForceRelativeDifference( const Vec3& f1, const Vec3& f2, double& forceNorm1, double& forceNorm2,
double& relativeDiff, FILE* log ) {
double diff = (f1[0] - f2[0])*(f1[0] - f2[0]) +
(f1[1] - f2[1])*(f1[1] - f2[1]) +
(f1[2] - f2[2])*(f1[2] - f2[2]);
forceNorm1 = sqrt( f1[0]*f1[0] + f1[1]*f1[1] + f1[2]*f1[2] );
forceNorm2 = sqrt( f2[0]*f2[0] + f2[1]*f2[1] + f2[2]*f2[2] );
if( forceNorm1 > 0.0 || forceNorm2 > 0.0 ){
relativeDiff = 2.0*sqrt( diff )/(forceNorm1+forceNorm2);
} else {
relativeDiff = 0.0;
}
return;
}
/**---------------------------------------------------------------------------------------
*
* Compare forces from two states
*
* @param state1 state1
* @param state2 state2
* @param relativeTolerance relative tolerance
* @param log if set, output forces
*
* @return number of entries with relative difference > tolerance
*
--------------------------------------------------------------------------------------- */
int compareForcesOfTwoStates( State& state1, State& state2, double relativeTolerance,
DoubleVector& stats, FILE* log ) {
int error = 0;
vector<Vec3> force1 = state1.getForces();
vector<Vec3> force2 = state2.getForces();
double averageRelativeDifference = 0.0;
double count = 0.0;
DoubleVector medians1( force1.size() );
DoubleVector medians2( force1.size() );
IntDoublePairVector relativeDifferences;
for( unsigned int ii = 0; ii < force1.size(); ii++ ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
getForceRelativeDifference( force1[ii], force2[ii], forceNorm1, forceNorm2, relativeDiff, log );
medians1[ii] = forceNorm1;
medians2[ii] = forceNorm2;
relativeDifferences.push_back( IntDoublePair(ii, relativeDiff ) );
averageRelativeDifference += relativeDiff;
count += 1.0;
if( relativeDiff > relativeTolerance ){
error++;
}
if( log ){
(void) fprintf( log, "F %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n", static_cast<unsigned int>(ii),
relativeDiff, force1[ii][0], force1[ii][1], force1[ii][2], force2[ii][0], force2[ii][1], force2[ii][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
// sort relative differences
std::sort( relativeDifferences.begin(), relativeDifferences.end(), TestIntDoublePair );
if( log ){
(void) fprintf( log, "\nEntries w/ largest relative differences.\n" );
for( unsigned int ii = relativeDifferences.size()-1; ii >= relativeDifferences.size()-10 && ii >= 0; ii-- ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
int index = relativeDifferences[ii].first;
getForceRelativeDifference( force1[index], force2[index], forceNorm1, forceNorm2, relativeDiff, log );
(void) fprintf( log, "Fs %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n",
static_cast<unsigned int>(index), relativeDiff,
force1[index][0], force1[index][1], force1[index][2],
force2[index][0], force2[index][1], force2[index][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
if( count > 0.0 ){
averageRelativeDifference /= count;
}
std::sort( medians1.begin(), medians1.end() );
std::sort( medians2.begin(), medians2.end() );
double median1 = medians1[medians1.size()/2];
double median2 = medians2[medians2.size()/2];
stats.resize( 4 );
stats[0] = averageRelativeDifference;
IntDoublePair pair = relativeDifferences[relativeDifferences.size()-1];
stats[1] = pair.second;
stats[2] = static_cast<double>(pair.first);
stats[3] = median1 < median2 ? median1 : median2;
return error;
}
void testEnergyEthaneSwitchingFunction( int useSwitchingFunction, FILE* log ) {
std::string methodName = "testEnergyEthaneSwitchingFunction";
System system;
const int numParticles = 8;
for( int i = 0; i < numParticles; i++ ){
system.addParticle(1.0);
}
double C_HBondDistance = 0.1097;
double C_CBondDistance = 0.1504;
NonbondedSoftcoreForce* nonbonded = new NonbondedSoftcoreForce();
nonbonded->setNonbondedMethod(NonbondedSoftcoreForce::NoCutoff);
double C_radius, C_gamma, C_charge, H_radius, H_gamma, H_charge;
int AM1_BCC = 1;
H_charge = -0.053;
C_charge = -3.0*H_charge;
if( AM1_BCC ){
C_radius = 0.180;
//C_radius = 0.360;
C_gamma = -0.2863;
C_gamma = 1.0;
H_radius = 0.125;
//H_radius = 0.25;
H_gamma = 0.2437;
H_gamma = 1.0;
//H_charge = C_charge = 0.0;
//H_gamma = C_gamma = 0.0;
} else {
C_radius = 0.215;
C_gamma = -1.1087;
H_radius = 0.150;
H_gamma = 0.1237;
}
// for ethane all Coulomb forces are excluded since all atoms 3 or
// fewer bonds away from all other atoms -- is this true for H's on
// difference carbons? -- should be computed in 14 ixn
double bornRadiusScaleFactorsEven = 0.5;
//double bornRadiusScaleFactorsEven = 1.0;
//double bornRadiusScaleFactorsOdd = 0.5;
double bornRadiusScaleFactorsOdd = 1.0;
if( log ){
(void) fprintf( log, "%s: Applying GB/VI\n", methodName.c_str() );
(void) fprintf( log, "C[%14.7e %14.7e %14.7e] H[%14.7e %14.7e %14.7e] scale[%.1f %.1f]\n",
C_charge, C_radius, C_gamma, H_charge, H_radius, H_gamma,
bornRadiusScaleFactorsEven, bornRadiusScaleFactorsOdd);
}
GBVISoftcoreForce* forceField = new GBVISoftcoreForce();
for( int i = 0; i < numParticles; i++ ){
forceField->addParticle( H_charge, H_radius, H_gamma, (i%2) ? bornRadiusScaleFactorsOdd : bornRadiusScaleFactorsEven);
nonbonded->addParticle( H_charge, H_radius, 0.0);
}
forceField->setParticleParameters( 1, C_charge, C_radius, C_gamma, bornRadiusScaleFactorsOdd);
nonbonded->setParticleParameters( 1, C_charge, C_radius, 0.0);
forceField->setParticleParameters( 4, C_charge, C_radius, C_gamma, bornRadiusScaleFactorsEven);
nonbonded->setParticleParameters( 4, C_charge, C_radius, 0.0);
// forceField->setParticleParameters( 8, C_charge, (C_radius+0.5), C_gamma, bornRadiusScaleFactorsEven);
// nonbonded->setParticleParameters( 8, C_charge, C_radius, 0.0);
if( useSwitchingFunction ){
forceField->setBornRadiusScalingMethod( GBVISoftcoreForce::QuinticSpline );
} else {
forceField->setBornRadiusScalingMethod( GBVISoftcoreForce::NoScaling );
}
forceField->addBond( 0, 1, C_HBondDistance );
forceField->addBond( 2, 1, C_HBondDistance );
forceField->addBond( 3, 1, C_HBondDistance );
forceField->addBond( 1, 4, C_CBondDistance );
forceField->addBond( 5, 4, C_HBondDistance );
forceField->addBond( 6, 4, C_HBondDistance );
forceField->addBond( 7, 4, C_HBondDistance );
std::vector<pair<int, int> > bonds;
std::vector<double> bondDistances;
bonds.push_back(pair<int, int>(0, 1));
bondDistances.push_back( C_HBondDistance );
bonds.push_back(pair<int, int>(2, 1));
bondDistances.push_back( C_HBondDistance );
bonds.push_back(pair<int, int>(3, 1));
bondDistances.push_back( C_HBondDistance );
bonds.push_back(pair<int, int>(1, 4));
bondDistances.push_back( C_CBondDistance );
bonds.push_back(pair<int, int>(5, 4));
bondDistances.push_back( C_HBondDistance );
bonds.push_back(pair<int, int>(6, 4));
bondDistances.push_back( C_HBondDistance );
bonds.push_back(pair<int, int>(7, 4));
bondDistances.push_back( C_HBondDistance );
nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
system.addForce(forceField);
system.addForce(nonbonded);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context referenceContext(system, integrator1, Platform::getPlatformByName( "Reference") );
Context context(system, integrator2, Platform::getPlatformByName( "Cuda") );
vector<Vec3> positions(numParticles);
positions[0] = Vec3(0.5480, 1.7661, 0.0000);
positions[1] = Vec3(0.7286, 0.8978, 0.6468);
positions[2] = Vec3(0.4974, 0.0000, 0.0588);
positions[3] = Vec3(0.0000, 0.9459, 1.4666);
positions[4] = Vec3(2.1421, 0.8746, 1.1615);
positions[5] = Vec3(2.3239, 0.0050, 1.8065);
positions[6] = Vec3(2.8705, 0.8295, 0.3416);
positions[7] = Vec3(2.3722, 1.7711, 1.7518);
//positions[8] = Vec3(2.1421, 0.8746, 2.1615);
vector<Vec3> originalPositions(numParticles);
for( int ii = 0; ii < numParticles; ii++ ){
originalPositions[ii][0] = positions[ii][0];
originalPositions[ii][1] = positions[ii][1];
originalPositions[ii][2] = positions[ii][2];
}
int tries = 1;
double positionIncrement = 0.15;
for( int ii = 0; ii < tries; ii++ ){
context.setPositions(positions);
referenceContext.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Forces | State::Energy);
if( log ){
(void) fprintf( log, "cudaE=%14.7e refE=%14.7e\n", state.getPotentialEnergy(), referenceState.getPotentialEnergy() );
}
// Take a small step in the direction of the energy gradient.
std::vector<double> stats;
if( compareForcesOfTwoStates( state, referenceState, 0.001, stats, log ) ){
ASSERT_EQUAL_TOL(0.0, 1.0, 0.01)
}
double norm = 0.0;
double forceSum[3] = { 0.0, 0.0, 0.0 };
for (int i = 0; i < numParticles; ++i) {
Vec3 f = state.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
forceSum[0] += f[0];
forceSum[1] += f[1];
forceSum[2] += f[2];
}
norm = std::sqrt(norm);
if( log ){
(void) fprintf( log, "Fsum [%14.7e %14.7e %14.7e] norm=%14.7e\n", forceSum[0], forceSum[1], forceSum[2], norm );
}
const double delta = 1e-03;
double step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
context.setPositions(positions);
State state2 = context.getState(State::Energy);
double diff = (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta;
double off = fabs( diff - norm )/norm;
if( log ){
(void) fprintf( log, "%2d Energies %.8e %.8e norms[%13.7e %13.7e] deltaNorms=%13.7e delta=%.2e\n",
ii, state.getPotentialEnergy(), state2.getPotentialEnergy(), diff, norm, off, delta );
}
// See whether the potential energy changed by the expected amount.
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta, 1e-3*abs(state.getPotentialEnergy()) );
if( ii < (tries-1) ){
for( int jj = 0; jj < numParticles; jj++ ){
positions[jj][0] = originalPositions[jj][0];
positions[jj][1] = originalPositions[jj][1];
positions[jj][2] = originalPositions[jj][2];
}
// positions[8][2] -= static_cast<double>(ii+1)*0.1;
// positions[8][2] -= 0.001;
if( log ){
(void) fprintf( log, "r48=%14.6e r28=%14.6e r24=%14.6e\n", positions[8][2]-positions[4][2], positions[8][2], positions[4][2] );
}
}
#if 0
int carbonIndex = 1;
int hydrogenIndex = 0;
while( hydrogenIndex < 8 ){
Vec3 carbonDelta;
for( int kk = 0; kk < 3; kk++ ){
positions[hydrogenIndex][kk] += positionIncrement*(positions[carbonIndex][kk] - positions[hydrogenIndex][kk] );
}
double dist = 0.0;
for( int kk = 0; kk < 3; kk++ ){
dist += (positions[carbonIndex][kk] - positions[hydrogenIndex][kk] )*(positions[carbonIndex][kk] - positions[hydrogenIndex][kk]);
}
(void) fprintf( log, "H=%d C=%d r=%14.6e\n", hydrogenIndex, carbonIndex, dist );
hydrogenIndex++;
if( hydrogenIndex == carbonIndex ){
hydrogenIndex++;
}
if( carbonIndex == 1 && hydrogenIndex == 4 ){
carbonIndex = 4;
hydrogenIndex = 5;
}
}
#endif
}
}
int main() {
try {
registerFreeEnergyCudaKernelFactories( );
//FILE* log = stderr;
FILE* log = NULL;
testSingleParticle( log );
testEnergyEthaneSwitchingFunction( 0, log );
testEnergyEthaneSwitchingFunction( 1, log );
} catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
plugins/freeEnergy/platforms/cuda/tests/TestCudaLJSoftcoreForce.cpp deleted 100644 → 0
View file @ c775bd19
/* -------------------------------------------------------------------------- *
* 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-2009 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 CustomGBForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "sfmt/SFMT.h"
#include "openmm/Context.h"
#include "openmm/CustomBondForce.h"
#include "openmm/CustomNonbondedForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/NonbondedSoftcoreForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include <iostream>
#include <vector>
#include <algorithm>
extern "C" void registerFreeEnergyCudaKernelFactories();
using namespace OpenMM;
using namespace std;
const double TOL = 1e-4;
static const int NoCutoff = 0;
static const int CutoffNonPeriodic = 1;
static const int CutoffPeriodic = 2;
void testNonbondedSoftcore( double lambda1, double lambda2, int nonbondedMethod, FILE* log ){
const int numMolecules = 70;
const int numParticles = numMolecules*2;
const double boxSize = 10.0;
const double reactionFieldDielectric = 80.0;
const double cutoffDistance = 0.4*boxSize;
// Create two systems: one with a NonbondedSoftcoreForce, and one using a CustomNonbondedForce to implement the same interaction.
System standardSystem;
System customSystem;
for (int i = 0; i < numParticles; i++) {
standardSystem.addParticle(1.0);
customSystem.addParticle(1.0);
}
standardSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
customSystem.setDefaultPeriodicBoxVectors( Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
NonbondedSoftcoreForce* nonbondedSoftcoreForce = new NonbondedSoftcoreForce();
CustomNonbondedForce* customNonbonded;
CustomBondForce* customBond;
if( nonbondedMethod == NoCutoff ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::NoCutoff );
customNonbonded = new CustomNonbondedForce("lambda*4*eps*(dem^2-dem)+138.935456*q/r;"
"q=q1*q2;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda);"
"sigma=0.5*(sigma1+sigma2);"
"eps=sqrt(eps1*eps2);"
"lambda=min(lambda1,lambda2)");
customNonbonded->setNonbondedMethod( CustomNonbondedForce::NoCutoff );
customBond = new CustomBondForce("lambda*4*eps*(dem^2-dem)+138.935456*q/r;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda)");
} else {
nonbondedSoftcoreForce->setCutoffDistance( cutoffDistance );
nonbondedSoftcoreForce->setReactionFieldDielectric( reactionFieldDielectric );
if( nonbondedMethod == CutoffNonPeriodic ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::CutoffNonPeriodic );
} else {
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::CutoffPeriodic );
}
customNonbonded = new CustomNonbondedForce("lambda*4*eps*(dem^2-dem)+138.935456*q*(1.0/r+(krf*r*r)-crf);"
"q=q1*q2;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda);"
"sigma=0.5*(sigma1+sigma2);"
"eps=sqrt(eps1*eps2);"
"lambda=min(lambda1,lambda2)");
customBond = new CustomBondForce("withinCutoff*(lambda*4*eps*(dem^2-dem)+138.935456*q*(1.0/r+(krf*r*r)-crf));"
"withinCutoff=step(cutoff-r);"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda)");
customNonbonded->setCutoffDistance( cutoffDistance );
if( nonbondedMethod == CutoffNonPeriodic ){
customNonbonded->setNonbondedMethod( CustomNonbondedForce::CutoffNonPeriodic );
} else {
customNonbonded->setNonbondedMethod( CustomNonbondedForce::CutoffPeriodic );
}
double eps2 = (reactionFieldDielectric - 1.0)/(2.0*reactionFieldDielectric+1.0);
double kValue = eps2/(cutoffDistance*cutoffDistance*cutoffDistance);
customNonbonded->addGlobalParameter("krf", kValue );
customBond->addGlobalParameter("krf", kValue );
double cValue = (1.0/cutoffDistance)*(3.0*reactionFieldDielectric)/(2.0*reactionFieldDielectric + 1.0);
customNonbonded->addGlobalParameter("crf", cValue );
customBond->addGlobalParameter("crf", cValue );
customBond->addGlobalParameter("cutoff", cutoffDistance );
}
customNonbonded->addPerParticleParameter("q");
customNonbonded->addPerParticleParameter("sigma");
customNonbonded->addPerParticleParameter("eps");
customNonbonded->addPerParticleParameter("lambda");
customBond->addPerBondParameter("q");
customBond->addPerBondParameter("sigma");
customBond->addPerBondParameter("eps");
customBond->addPerBondParameter("lambda");
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<double> params(4);
// periodic boundary conditions not possible w/ CustomBond?
int includeExceptions = nonbondedMethod != NoCutoff ? 0 : 1;
for (int i = 0; i < numMolecules; i++) {
if (i < numMolecules/2) {
double charge = 1.0;
double sigma1 = 0.2;
double sigma2 = 0.2;
double eps1 = 0.5;
double eps2 = 0.5;
//eps1 = eps2 = 0.0;
nonbondedSoftcoreForce->addParticle( charge, sigma1, eps1, lambda1);
nonbondedSoftcoreForce->addParticle(-charge, sigma2, eps2, lambda1);
params[0] = charge;
params[1] = sigma1;
params[2] = eps1;
params[3] = lambda1;
customNonbonded->addParticle(params);
params[0] = -charge;
params[1] = sigma2;
params[2] = eps2;
customNonbonded->addParticle(params);
if( includeExceptions && i && ((i%4) == 0) ){
vector<double> bondParams(4);
nonbondedSoftcoreForce->addException(i-4, i, charge*charge, sigma1, eps1, false, lambda1);
customNonbonded->addExclusion( i-4,i);
bondParams[0] = charge*charge;
bondParams[1] = sigma1;
bondParams[2] = eps1;
bondParams[3] = lambda1;
customBond->addBond(i-4,i, bondParams );
}
} else {
double charge = 1.0;
double sigma1 = 0.2;
double sigma2 = 0.1;
double eps1 = 0.8;
double eps2 = 0.8;
//eps1 = eps2 = 0.0;
nonbondedSoftcoreForce->addParticle( charge, sigma1, eps1, lambda2);
nonbondedSoftcoreForce->addParticle(-charge, sigma2, eps2, lambda2);
params[0] = charge;
params[1] = sigma1;
params[2] = eps1;
params[3] = lambda2;
customNonbonded->addParticle(params);
params[0] = -charge;
params[1] = sigma2;
params[2] = eps2;
customNonbonded->addParticle(params);
if( includeExceptions && i && ((i%4) == 0) ){
vector<double> bondParams(4);
nonbondedSoftcoreForce->addException(i-4, i, charge*charge, sigma1, eps1, false, lambda2);
customNonbonded->addExclusion( i-4,i);
bondParams[0] = charge*charge;
bondParams[1] = sigma1;
bondParams[2] = eps1;
bondParams[3] = lambda2;
customBond->addBond(i-4,i, bondParams );
}
}
positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
velocities[2*i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
velocities[2*i+1] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
}
standardSystem.addForce(nonbondedSoftcoreForce);
customSystem.addForce(customNonbonded);
customSystem.addForce(customBond);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context context1(standardSystem, integrator1, Platform::getPlatformByName( "Cuda"));
context1.setPositions(positions);
context1.setVelocities(velocities);
State state1 = context1.getState(State::Forces | State::Energy);
Context context2(customSystem, integrator2, Platform::getPlatformByName( "Reference"));
context2.setPositions(positions);
context2.setVelocities(velocities);
State state2 = context2.getState(State::Forces | State::Energy);
double diff = 0.0;
static int previousNonbondedMethod = -1;
if( state1.getPotentialEnergy() - state2.getPotentialEnergy() != 0.0 ){
diff = fabs( state1.getPotentialEnergy() - state2.getPotentialEnergy() )/( fabs( state1.getPotentialEnergy() ) + fabs( state2.getPotentialEnergy() ) );
}
if( previousNonbondedMethod != nonbondedMethod ){
if( log ){
(void) fprintf( log, "\n" );
}
previousNonbondedMethod = nonbondedMethod;
}
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
double maxDiff = -1.0;
int maxDiffIndex = -1;
for (int i = 0; i < numParticles; i++) {
Vec3 f1 = state1.getForces()[i];
Vec3 f2 = state2.getForces()[i];
double f1N = sqrt( (f1[0]*f1[0]) + (f1[1]*f1[1]) + (f1[2]*f1[2]) );
double f2N = sqrt( (f2[0]*f2[0]) + (f2[1]*f2[1]) + (f2[2]*f2[2]) );
double diff = (f1[0]-f2[0])*(f1[0]-f2[0]) +
(f1[1]-f2[1])*(f1[1]-f2[1]) +
(f1[2]-f2[2])*(f1[2]-f2[2]);
if( f1N > 0.0 || f1N > 0.0 ){
diff = 2.0*sqrt( diff )/(f1N + f2N);
}
if( diff > maxDiff ){
maxDiff = diff;
maxDiffIndex = i;
}
/*
(void) fprintf( log, "%4d %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] \n", i, diff,
f1[0], f1[1],f1[2], f2[0], f2[1],f2[2] );
*/
// ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
ASSERT( diff < 1.0e-3);
}
if( log ){
(void) fprintf( log, "%d %d %10.1f %10.1f %15.7e %15.7e %15.7e maxFDiff=%15.7e %d\n",
nonbondedMethod, includeExceptions, lambda1, lambda2, diff,
state1.getPotentialEnergy(), state2.getPotentialEnergy(), maxDiff, maxDiffIndex); fflush( log );
}
}
int main() {
try {
registerFreeEnergyCudaKernelFactories( );
// test various combinations of lambdas and boundary conditions/cutoffs
FILE* log = NULL;
testNonbondedSoftcore( 1.0, 1.0 , NoCutoff, log );
testNonbondedSoftcore( 1.0, 0.0 , NoCutoff, log );
testNonbondedSoftcore( 1.0, 0.5 , NoCutoff, log );
testNonbondedSoftcore( 0.0, 0.0 , NoCutoff, log );
testNonbondedSoftcore( 1.0, 1.0 , CutoffNonPeriodic, log );
testNonbondedSoftcore( 1.0, 0.0 , CutoffNonPeriodic, log );
testNonbondedSoftcore( 1.0, 0.5 , CutoffNonPeriodic, log );
testNonbondedSoftcore( 0.0, 0.0 , CutoffNonPeriodic, log );
testNonbondedSoftcore( 1.0, 1.0 , CutoffPeriodic, log );
testNonbondedSoftcore( 1.0, 0.0 , CutoffPeriodic, log );
testNonbondedSoftcore( 1.0, 0.5 , CutoffPeriodic, log );
testNonbondedSoftcore( 0.0, 0.0 , CutoffPeriodic, log );
} catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
plugins/freeEnergy/platforms/cuda/tests/TestCudaOBCSoftcoreForce.cpp deleted 100644 → 0
View file @ c775bd19
/* -------------------------------------------------------------------------- *
* 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-2009 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 implementations of GBSAOBCSoftcoreForce
*/
#include "openmm/System.h"
#include "openmm/internal/AssertionUtilities.h"
#include "sfmt/SFMT.h"
#include "openmm/Context.h"
#include "openmm/NonbondedSoftcoreForce.h"
#include "openmm/GBSAOBCSoftcoreForce.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "openmm/VerletIntegrator.h"
#include <iostream>
#include <vector>
#include <algorithm>
#include <iomanip>
typedef std::vector<double> DoubleVector;
typedef DoubleVector::iterator DoubleVectorI;
typedef DoubleVector::const_iterator DoubleVectorCI;
typedef std::vector<DoubleVector> DoubleVectorVector;
typedef std::pair<int, double> IntDoublePair;
typedef std::vector<IntDoublePair> IntDoublePairVector;
typedef IntDoublePairVector::iterator IntDoublePairVectorI;
typedef IntDoublePairVector::const_iterator IntDoublePairVectorCI;
extern "C" void registerFreeEnergyCudaKernelFactories();
using namespace OpenMM;
using namespace std;
void testSingleParticle( FILE* log ) {
System system;
system.addParticle(2.0);
VerletIntegrator integrator(0.01);
GBSAOBCSoftcoreForce* forceField = new GBSAOBCSoftcoreForce();
double charge = 0.5;
double radius = 0.15;
double scaleFactor = 1.0;
forceField->addParticle(charge, radius, scaleFactor);
system.addForce(forceField);
NonbondedSoftcoreForce* nonbonded = new NonbondedSoftcoreForce();
nonbonded->setNonbondedMethod(NonbondedSoftcoreForce::NoCutoff);
nonbonded->addParticle( charge, 1.0, 0.0);
system.addForce(nonbonded);
Context context(system, integrator, Platform::getPlatformByName( "Cuda") );
vector<Vec3> positions(1);
positions[0] = Vec3(0, 0, 0);
context.setPositions(positions);
State state = context.getState(State::Energy);
double bornRadius = radius-0.009;
double eps0 = EPSILON0;
double tau = (1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric());
double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric())/bornRadius;
double extendedRadius = bornRadius+0.14; // probe radius
double nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.15/bornRadius, 6.0); // Where did this formula come from? Just copied it from CpuImplicitSolvent.cpp
double expectedE = (bornEnergy+nonpolarEnergy);
double obtainedE = state.getPotentialEnergy();
double diff = fabs( obtainedE - expectedE );
if( log ){
(void) fprintf( log, "testSingleParticle expected=%14.6e obtained=%14.6e diff=%14.6e breakdown:[%14.6e %14.6e]\n",
expectedE, obtainedE, diff, bornEnergy, nonpolarEnergy );
}
ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
}
/**
* Predicate for sorting <int,double> pair
*
* @param d1 first IntDoublePair to compare
* @param d2 second IntDoublePair to compare
*
*/
bool TestIntDoublePair( const IntDoublePair& d1, const IntDoublePair& d2 ){
return d1.second < d2.second;
}
/**---------------------------------------------------------------------------------------
*
* Get relative difference between two forces
*
* @param f1 force1
* @param f2 force2
* @param forceNorm1 output norm of force1
* @param forceNorm2 output norm of force2
* @param relativeDiff output relative difference between force norms
* @param log if set, output forces
*
*
--------------------------------------------------------------------------------------- */
static void getForceRelativeDifference( const Vec3& f1, const Vec3& f2, double& forceNorm1, double& forceNorm2,
double& relativeDiff, FILE* log ) {
double diff = (f1[0] - f2[0])*(f1[0] - f2[0]) +
(f1[1] - f2[1])*(f1[1] - f2[1]) +
(f1[2] - f2[2])*(f1[2] - f2[2]);
forceNorm1 = sqrt( f1[0]*f1[0] + f1[1]*f1[1] + f1[2]*f1[2] );
forceNorm2 = sqrt( f2[0]*f2[0] + f2[1]*f2[1] + f2[2]*f2[2] );
if( forceNorm1 > 0.0 || forceNorm2 > 0.0 ){
relativeDiff = 2.0*sqrt( diff )/(forceNorm1+forceNorm2);
} else {
relativeDiff = 0.0;
}
return;
}
/**---------------------------------------------------------------------------------------
*
* Compare forces from two states
*
* @param state1 state1
* @param state2 state2
* @param relativeTolerance relative tolerance
* @param log if set, output forces
*
* @return number of entries with relative difference > tolerance
*
--------------------------------------------------------------------------------------- */
int compareForcesOfTwoStates( State& state1, State& state2, double relativeTolerance,
DoubleVector& stats, FILE* log ) {
int error = 0;
vector<Vec3> force1 = state1.getForces();
vector<Vec3> force2 = state2.getForces();
double averageRelativeDifference = 0.0;
double count = 0.0;
DoubleVector medians1( force1.size() );
DoubleVector medians2( force1.size() );
IntDoublePairVector relativeDifferences;
for( unsigned int ii = 0; ii < force1.size(); ii++ ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
getForceRelativeDifference( force1[ii], force2[ii], forceNorm1, forceNorm2, relativeDiff, log );
medians1[ii] = forceNorm1;
medians2[ii] = forceNorm2;
relativeDifferences.push_back( IntDoublePair(ii, relativeDiff ) );
averageRelativeDifference += relativeDiff;
count += 1.0;
if( relativeDiff > relativeTolerance ){
error++;
}
if( log ){
(void) fprintf( log, "F %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n", static_cast<unsigned int>(ii),
relativeDiff, force1[ii][0], force1[ii][1], force1[ii][2], force2[ii][0], force2[ii][1], force2[ii][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
// sort relative differences
std::sort( relativeDifferences.begin(), relativeDifferences.end(), TestIntDoublePair );
if( log ){
(void) fprintf( log, "\nEntries w/ largest relative differences.\n" );
for( unsigned int ii = relativeDifferences.size()-1; ii >= relativeDifferences.size()-10 && ii >= 0; ii-- ){
double forceNorm1;
double forceNorm2;
double relativeDiff;
int index = relativeDifferences[ii].first;
getForceRelativeDifference( force1[index], force2[index], forceNorm1, forceNorm2, relativeDiff, log );
(void) fprintf( log, "Fs %6u %15.7e [%15.7e %15.7e %15.7e] [%15.7e %15.7e %15.7e] %15.7e %15.7e %s\n",
static_cast<unsigned int>(index), relativeDiff,
force1[index][0], force1[index][1], force1[index][2],
force2[index][0], force2[index][1], force2[index][2],
forceNorm1, forceNorm2, (relativeDiff < relativeTolerance ? "":"XXXXXX") );
}
}
if( count > 0.0 ){
averageRelativeDifference /= count;
}
std::sort( medians1.begin(), medians1.end() );
std::sort( medians2.begin(), medians2.end() );
double median1 = medians1[medians1.size()/2];
double median2 = medians2[medians2.size()/2];
stats.resize( 4 );
stats[0] = averageRelativeDifference;
IntDoublePair pair = relativeDifferences[relativeDifferences.size()-1];
stats[1] = pair.second;
stats[2] = static_cast<double>(pair.first);
stats[3] = median1 < median2 ? median1 : median2;
return error;
}
void testEnergyEthaneSwitchingFunction( int useSwitchingFunction, FILE* log ) {
std::string methodName = "testEnergyEthaneSwitchingFunction";
System system;
const int numParticles = 8;
for( int i = 0; i < numParticles; i++ ){
system.addParticle(1.0);
}
double C_HBondDistance = 0.1097;
double C_CBondDistance = 0.1504;
NonbondedSoftcoreForce* nonbonded = new NonbondedSoftcoreForce();
nonbonded->setNonbondedMethod(NonbondedSoftcoreForce::NoCutoff);
double C_radius, C_gamma, C_charge, H_radius, H_gamma, H_charge;
int AM1_BCC = 1;
H_charge = -0.053;
C_charge = -3.0*H_charge;
if( AM1_BCC ){
C_radius = 0.180;
//C_radius = 0.360;
C_gamma = -0.2863;
C_gamma = 1.0;
H_radius = 0.125;
//H_radius = 0.25;
H_gamma = 0.2437;
H_gamma = 1.0;
//H_charge = C_charge = 0.0;
//H_gamma = C_gamma = 0.0;
} else {
C_radius = 0.215;
C_gamma = -1.1087;
H_radius = 0.150;
H_gamma = 0.1237;
}
// for ethane all Coulomb forces are excluded since all atoms 3 or
// fewer bonds away from all other atoms -- is this true for H's on
// difference carbons? -- should be computed in 14 ixn
double bornRadiusScaleFactorsEven = 0.5;
//double bornRadiusScaleFactorsEven = 1.0;
//double bornRadiusScaleFactorsOdd = 0.5;
double bornRadiusScaleFactorsOdd = 1.0;
if( log ){
(void) fprintf( log, "%s: Applying GB/VI\n", methodName.c_str() );
(void) fprintf( log, "C[%14.7e %14.7e %14.7e] H[%14.7e %14.7e %14.7e] scale[%.1f %.1f]\n",
C_charge, C_radius, C_gamma, H_charge, H_radius, H_gamma,
bornRadiusScaleFactorsEven, bornRadiusScaleFactorsOdd);
}
GBSAOBCSoftcoreForce* forceField = new GBSAOBCSoftcoreForce();
for( int i = 0; i < numParticles; i++ ){
forceField->addParticle( H_charge, H_radius, H_gamma, (i%2) ? bornRadiusScaleFactorsOdd : bornRadiusScaleFactorsEven);
nonbonded->addParticle( H_charge, H_radius, 0.0);
}
forceField->setParticleParameters( 1, C_charge, C_radius, C_gamma, bornRadiusScaleFactorsOdd);
nonbonded->setParticleParameters( 1, C_charge, C_radius, 0.0);
forceField->setParticleParameters( 4, C_charge, C_radius, C_gamma, bornRadiusScaleFactorsEven);
nonbonded->setParticleParameters( 4, C_charge, C_radius, 0.0);
// forceField->setParticleParameters( 8, C_charge, (C_radius+0.5), C_gamma, bornRadiusScaleFactorsEven);
// nonbonded->setParticleParameters( 8, C_charge, C_radius, 0.0);
std::vector<pair<int, int> > bonds;
std::vector<double> bondDistances;
bonds.push_back(pair<int, int>(0, 1));
bondDistances.push_back( C_HBondDistance );
bonds.push_back(pair<int, int>(2, 1));
bondDistances.push_back( C_HBondDistance );
bonds.push_back(pair<int, int>(3, 1));
bondDistances.push_back( C_HBondDistance );
bonds.push_back(pair<int, int>(1, 4));
bondDistances.push_back( C_CBondDistance );
bonds.push_back(pair<int, int>(5, 4));
bondDistances.push_back( C_HBondDistance );
bonds.push_back(pair<int, int>(6, 4));
bondDistances.push_back( C_HBondDistance );
bonds.push_back(pair<int, int>(7, 4));
bondDistances.push_back( C_HBondDistance );
nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
system.addForce(forceField);
system.addForce(nonbonded);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context referenceContext(system, integrator1, Platform::getPlatformByName( "Reference") );
Context context(system, integrator2, Platform::getPlatformByName( "Cuda") );
vector<Vec3> positions(numParticles);
positions[0] = Vec3(0.5480, 1.7661, 0.0000);
positions[1] = Vec3(0.7286, 0.8978, 0.6468);
positions[2] = Vec3(0.4974, 0.0000, 0.0588);
positions[3] = Vec3(0.0000, 0.9459, 1.4666);
positions[4] = Vec3(2.1421, 0.8746, 1.1615);
positions[5] = Vec3(2.3239, 0.0050, 1.8065);
positions[6] = Vec3(2.8705, 0.8295, 0.3416);
positions[7] = Vec3(2.3722, 1.7711, 1.7518);
//positions[8] = Vec3(2.1421, 0.8746, 2.1615);
vector<Vec3> originalPositions(numParticles);
for( int ii = 0; ii < numParticles; ii++ ){
originalPositions[ii][0] = positions[ii][0];
originalPositions[ii][1] = positions[ii][1];
originalPositions[ii][2] = positions[ii][2];
}
int tries = 1;
double positionIncrement = 0.15;
for( int ii = 0; ii < tries; ii++ ){
context.setPositions(positions);
referenceContext.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Forces | State::Energy);
if( log ){
(void) fprintf( log, "cudaE=%14.7e refE=%14.7e\n", state.getPotentialEnergy(), referenceState.getPotentialEnergy() );
}
// Take a small step in the direction of the energy gradient.
DoubleVector stats;
if( compareForcesOfTwoStates( state, referenceState, 0.001, stats, log ) ){
ASSERT_EQUAL_TOL(0.0, 1.0, 0.01)
}
double norm = 0.0;
double forceSum[3] = { 0.0, 0.0, 0.0 };
for (int i = 0; i < numParticles; ++i) {
Vec3 f = state.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
forceSum[0] += f[0];
forceSum[1] += f[1];
forceSum[2] += f[2];
}
norm = std::sqrt(norm);
if( log ){
(void) fprintf( log, "Fsum [%14.7e %14.7e %14.7e] norm=%14.7e\n", forceSum[0], forceSum[1], forceSum[2], norm );
}
const double delta = 1e-03;
double step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
context.setPositions(positions);
State state2 = context.getState(State::Energy);
double diff = (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta;
double off = fabs( diff - norm )/norm;
if( log ){
(void) fprintf( log, "%2d Energies %.8e %.8e norms[%13.7e %13.7e] deltaNorms=%13.7e delta=%.2e\n",
ii, state.getPotentialEnergy(), state2.getPotentialEnergy(), diff, norm, off, delta );
}
// See whether the potential energy changed by the expected amount.
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state.getPotentialEnergy())/delta, 1e-3*abs(state.getPotentialEnergy()) );
if( ii < (tries-1) ){
for( int jj = 0; jj < numParticles; jj++ ){
positions[jj][0] = originalPositions[jj][0];
positions[jj][1] = originalPositions[jj][1];
positions[jj][2] = originalPositions[jj][2];
}
// positions[8][2] -= static_cast<double>(ii+1)*0.1;
// positions[8][2] -= 0.001;
if( log ){
(void) fprintf( log, "r48=%14.6e r28=%14.6e r24=%14.6e\n", positions[8][2]-positions[4][2], positions[8][2], positions[4][2] );
}
}
}
}
int main() {
try {
registerFreeEnergyCudaKernelFactories( );
//FILE* log = stderr;
FILE* log = NULL;
testSingleParticle( log );
testEnergyEthaneSwitchingFunction( 0, log );
testEnergyEthaneSwitchingFunction( 1, log );
} catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
plugins/freeEnergy/platforms/reference/include/ReferenceFreeEnergyKernelFactory.h deleted 100644 → 0
View file @ c775bd19
#ifndef OPENMM_REFERENCE_FREE_ENERGY_KERNEL_FACTORY_H_
#define OPENMM_REFERENCE_FREE_ENERGY_KERNEL_FACTORY_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/KernelFactory.h"
namespace OpenMM {
/**
* This KernelFactory creates all kernels for ReferencePlatform.
*/
class OPENMM_EXPORT ReferenceFreeEnergyKernelFactory : public KernelFactory {
public:
KernelImpl* createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const;
};
} // namespace OpenMM
#endif /*OPENMM_REFERENCE_FREE_ENERGY_KERNEL_FACTORY_H_*/
plugins/freeEnergy/platforms/reference/src/ReferenceFreeEnergyKernelFactory.cpp deleted 100644 → 0
View file @ c775bd19
/* -------------------------------------------------------------------------- *
* 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. *
* -------------------------------------------------------------------------- */
#include "ReferenceFreeEnergyKernelFactory.h"
#include "ReferenceFreeEnergyKernels.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/OpenMMException.h"
using namespace OpenMM;
#if defined(WIN32)
#include <windows.h>
extern "C" void initFreeEnergyReferenceKernels();
BOOL WINAPI DllMain(HANDLE hModule, DWORD ul_reason_for_call, LPVOID lpReserved) {
if (ul_reason_for_call == DLL_PROCESS_ATTACH)
initFreeEnergyReferenceKernels();
return TRUE;
}
#else
extern "C" void __attribute__((constructor)) initFreeEnergyReferenceKernels();
#endif
extern "C" void initFreeEnergyReferenceKernels() {
for( int ii = 0; ii < Platform::getNumPlatforms(); ii++ ){
Platform& platform = Platform::getPlatform(ii);
if( platform.getName().compare( "Reference" ) == 0 ){
ReferenceFreeEnergyKernelFactory* factory = new ReferenceFreeEnergyKernelFactory();
platform.registerKernelFactory(CalcNonbondedSoftcoreForceKernel::Name(), factory);
platform.registerKernelFactory(CalcGBSAOBCSoftcoreForceKernel::Name(), factory);
platform.registerKernelFactory(CalcGBVISoftcoreForceKernel::Name(), factory);
}
}
}
KernelImpl* ReferenceFreeEnergyKernelFactory::createKernelImpl(std::string name, const Platform& platform, ContextImpl& context) const {
ReferencePlatform::PlatformData& data = *static_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
if (name == CalcNonbondedSoftcoreForceKernel::Name())
return new ReferenceFreeEnergyCalcNonbondedSoftcoreForceKernel(name, platform);
if (name == CalcGBSAOBCSoftcoreForceKernel::Name())
return new ReferenceFreeEnergyCalcGBSAOBCSoftcoreForceKernel(name, platform);
if (name == CalcGBVISoftcoreForceKernel::Name())
return new ReferenceFreeEnergyCalcGBVISoftcoreForceKernel(name, platform);
throw OpenMMException( (std::string("Tried to create kernel with illegal kernel name '") + name + "'").c_str() );
}
plugins/freeEnergy/platforms/reference/src/ReferenceFreeEnergyKernels.cpp deleted 100644 → 0
View file @ c775bd19
/* -------------------------------------------------------------------------- *
* 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-2009 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ReferenceFreeEnergyKernels.h"
#include "gbsa/CpuGBVISoftcore.h"
#include "gbsa/CpuObcSoftcore.h"
#include "SimTKReference/ReferenceFreeEnergyLJCoulomb14Softcore.h"
#include "SimTKReference/ReferenceFreeEnergyLJCoulombSoftcoreIxn.h"
#include "ReferenceBondForce.h"
#include "openmm/System.h"
#include "openmm/internal/ContextImpl.h"
#include "SimTKUtilities/SimTKOpenMMUtilities.h"
#include <cmath>
#include <limits>
using namespace std;
using namespace OpenMM;
static int** allocateIntArray(int length, int width) {
int** array = new int*[length];
for (int i = 0; i < length; ++i)
array[i] = new int[width];
return array;
}
static RealOpenMM** allocateRealArray(int length, int width) {
RealOpenMM** array = new RealOpenMM*[length];
for (int i = 0; i < length; ++i)
array[i] = new RealOpenMM[width];
return array;
}
static void disposeIntArray(int** array, int size) {
if (array) {
for (int i = 0; i < size; ++i)
delete[] array[i];
delete[] array;
}
}
static void disposeRealArray(RealOpenMM** array, int size) {
if (array) {
for (int i = 0; i < size; ++i)
delete[] array[i];
delete[] array;
}
}
static vector<RealVec>& extractPositions(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<RealVec>*) data->positions);
}
static vector<RealVec>& extractVelocities(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<RealVec>*) data->velocities);
}
static vector<RealVec>& extractForces(ContextImpl& context) {
ReferencePlatform::PlatformData* data = reinterpret_cast<ReferencePlatform::PlatformData*>(context.getPlatformData());
return *((vector<RealVec>*) data->forces);
}
ReferenceFreeEnergyCalcNonbondedSoftcoreForceKernel::~ReferenceFreeEnergyCalcNonbondedSoftcoreForceKernel() {
disposeRealArray(particleParamArray, numParticles);
disposeIntArray(exclusionArray, numParticles);
disposeIntArray(bonded14IndexArray, num14);
disposeRealArray(bonded14ParamArray, num14);
if (neighborList != NULL)
delete neighborList;
}
void ReferenceFreeEnergyCalcNonbondedSoftcoreForceKernel::initialize(const System& system, const NonbondedSoftcoreForce& force) {
// Identify which exceptions are 1-4 interactions.
numParticles = force.getNumParticles();
exclusions.resize(numParticles);
std::vector<int> nb14s;
for (int i = 0; i < force.getNumExceptions(); i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon, softcoreLJLambda;
force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon, softcoreLJLambda);
exclusions[particle1].insert(particle2);
exclusions[particle2].insert(particle1);
if (chargeProd != 0.0 || epsilon != 0.0)
nb14s.push_back(i);
}
// Build the arrays.
particleParamArray = allocateRealArray(numParticles, 4);
RealOpenMM sqrtEps = static_cast<RealOpenMM>( std::sqrt(138.935485) );
for (int i = 0; i < numParticles; ++i) {
double charge, radius, depth, softcoreLJLambda;
force.getParticleParameters(i, charge, radius, depth,softcoreLJLambda);
particleParamArray[i][0] = static_cast<RealOpenMM>(0.5*radius);
particleParamArray[i][1] = static_cast<RealOpenMM>(2.0*sqrt(depth));
particleParamArray[i][2] = static_cast<RealOpenMM>(charge*sqrtEps);
particleParamArray[i][3] = static_cast<RealOpenMM>(softcoreLJLambda);
}
this->exclusions = exclusions;
exclusionArray = new int*[numParticles];
for (int i = 0; i < numParticles; ++i) {
exclusionArray[i] = new int[exclusions[i].size()+1];
exclusionArray[i][0] = exclusions[i].size();
int index = 0;
for (std::set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
exclusionArray[i][++index] = *iter;
}
num14 = nb14s.size();
bonded14IndexArray = allocateIntArray(num14, 2);
bonded14ParamArray = allocateRealArray(num14, 4);
for (int i = 0; i < num14; ++i) {
int particle1, particle2;
double charge, radius, depth, softcoreLJLambda;
force.getExceptionParameters(nb14s[i], particle1, particle2, charge, radius, depth, softcoreLJLambda);
bonded14IndexArray[i][0] = particle1;
bonded14IndexArray[i][1] = particle2;
bonded14ParamArray[i][0] = static_cast<RealOpenMM>(radius);
bonded14ParamArray[i][1] = static_cast<RealOpenMM>(4.0*depth);
bonded14ParamArray[i][2] = static_cast<RealOpenMM>(charge*sqrtEps*sqrtEps);
bonded14ParamArray[i][3] = static_cast<RealOpenMM>(softcoreLJLambda);
}
nonbondedMethod = CalcNonbondedSoftcoreForceKernel::NonbondedSoftcoreMethod(force.getNonbondedMethod());
nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
Vec3 boxVectors[3];
system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
if (nonbondedMethod == NoCutoff)
neighborList = NULL;
else
neighborList = new NeighborList();
rfDielectric = (RealOpenMM)force.getReactionFieldDielectric();
}
double ReferenceFreeEnergyCalcNonbondedSoftcoreForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
RealOpenMM energy = 0;
ReferenceFreeEnergyLJCoulombSoftcoreIxn clj;
bool periodic = (nonbondedMethod == CutoffPeriodic);
if (nonbondedMethod != NoCutoff) {
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodicBoxSize, periodic, nonbondedCutoff, 0.0);
clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
}
if (periodic)
clj.setPeriodic(periodicBoxSize);
clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
ReferenceBondForce refBondForce;
ReferenceFreeEnergyLJCoulomb14Softcore nonbonded14;
if (nonbondedMethod == CutoffNonPeriodic || nonbondedMethod == CutoffPeriodic)
nonbonded14.setUseCutoff(nonbondedCutoff, rfDielectric);
refBondForce.calculateForce(num14, bonded14IndexArray, posData, bonded14ParamArray, forceData, &energy, nonbonded14);
return energy;
}
ReferenceFreeEnergyCalcGBSAOBCSoftcoreForceKernel::~ReferenceFreeEnergyCalcGBSAOBCSoftcoreForceKernel() {
if (obc) {
ObcSoftcoreParameters* obcSoftcoreParameters = obc->getObcSoftcoreParameters();
delete obcSoftcoreParameters;
delete obc;
}
}
void ReferenceFreeEnergyCalcGBSAOBCSoftcoreForceKernel::initialize(const System& system, const GBSAOBCSoftcoreForce& force) {
int numParticles = system.getNumParticles();
charges.resize(numParticles);
std::vector<RealOpenMM> atomicRadii(numParticles);
std::vector<RealOpenMM> scaleFactors(numParticles);
std::vector<RealOpenMM> nonPolarScaleFactors(numParticles);
for (int i = 0; i < numParticles; ++i) {
double charge, radius, scalingFactor, nonPolarScaleFactor;
force.getParticleParameters(i, charge, radius, scalingFactor, nonPolarScaleFactor);
charges[i] = static_cast<RealOpenMM>(charge);
atomicRadii[i] = static_cast<RealOpenMM>(radius);
scaleFactors[i] = static_cast<RealOpenMM>(scalingFactor);
nonPolarScaleFactors[i] = static_cast<RealOpenMM>(nonPolarScaleFactor);
}
ObcSoftcoreParameters* obcParameters = new ObcSoftcoreParameters(numParticles, ObcSoftcoreParameters::ObcTypeII);
obcParameters->setAtomicRadii(atomicRadii);
obcParameters->setScaledRadiusFactors(scaleFactors);
obcParameters->setNonPolarScaleFactors(nonPolarScaleFactors);
obcParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
obcParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
// nonPolarPrefactor is in units of kJ/mol/nm^2 convert to kcal/mol/A^2
// to be consistent w/ polar part of OBC calculation
RealOpenMM nonPolarPrefactor = static_cast<RealOpenMM>(force.getNonPolarPrefactor());
nonPolarPrefactor /= static_cast<RealOpenMM>(418.4);
obcParameters->setNonPolarPrefactor( nonPolarPrefactor );
// If there is a NonbondedForce in this system, use it to initialize cutoffs and periodic boundary conditions.
for (int i = 0; i < system.getNumForces(); i++) {
const NonbondedSoftcoreForce* nonbonded = dynamic_cast<const NonbondedSoftcoreForce*>(&system.getForce(i));
if (nonbonded != NULL) {
if (nonbonded->getNonbondedMethod() != NonbondedSoftcoreForce::NoCutoff){
obcParameters->setUseCutoff(static_cast<RealOpenMM>(nonbonded->getCutoffDistance()));
}
if (nonbonded->getNonbondedMethod() == NonbondedSoftcoreForce::CutoffPeriodic) {
Vec3 boxVectors[3];
system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
RealOpenMM periodicBoxSize[3];
periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
obcParameters->setPeriodic(periodicBoxSize);
}
break;
}
}
obc = new CpuObcSoftcore(obcParameters);
obc->setIncludeAceApproximation(true);
}
double ReferenceFreeEnergyCalcGBSAOBCSoftcoreForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
return obc->computeBornEnergyForces(posData, charges, forceData);
}
ReferenceFreeEnergyCalcGBVISoftcoreForceKernel::~ReferenceFreeEnergyCalcGBVISoftcoreForceKernel() {
if (gbviSoftcore) {
delete gbviSoftcore;
}
}
void ReferenceFreeEnergyCalcGBVISoftcoreForceKernel::initialize(const System& system, const GBVISoftcoreForce& force, const std::vector<double> & inputScaledRadii ) {
int numParticles = system.getNumParticles();
charges.resize(numParticles);
std::vector<RealOpenMM> atomicRadii(numParticles);
std::vector<RealOpenMM> scaledRadii(numParticles);
std::vector<RealOpenMM> gammas(numParticles);
std::vector<RealOpenMM> bornRadiusScaleFactors(numParticles);
for (int i = 0; i < numParticles; ++i) {
double charge, radius, gamma, bornRadiusScaleFactor;
force.getParticleParameters(i, charge, radius, gamma, bornRadiusScaleFactor);
charges[i] = static_cast<RealOpenMM>(charge);
atomicRadii[i] = static_cast<RealOpenMM>(radius);
gammas[i] = static_cast<RealOpenMM>(gamma);
scaledRadii[i] = static_cast<RealOpenMM>(inputScaledRadii[i]);
bornRadiusScaleFactors[i] = static_cast<RealOpenMM>(bornRadiusScaleFactor);
}
GBVISoftcoreParameters* gBVIParameters = new GBVISoftcoreParameters(numParticles);
gBVIParameters->setAtomicRadii(atomicRadii);
gBVIParameters->setGammaParameters(gammas);
gBVIParameters->setBornRadiusScaleFactors(bornRadiusScaleFactors);
gBVIParameters->setScaledRadii(scaledRadii);
// switching function/scaling
// quintic spline
if( force.getBornRadiusScalingMethod() == GBVISoftcoreForce::QuinticSpline ){
gBVIParameters->setBornRadiusScalingSoftcoreMethod( GBVISoftcoreParameters::QuinticSpline );
gBVIParameters->setQuinticLowerLimitFactor( static_cast<RealOpenMM>(force.getQuinticLowerLimitFactor()) );
gBVIParameters->setQuinticUpperBornRadiusLimit( static_cast<RealOpenMM>(force.getQuinticUpperBornRadiusLimit()) );
}
gBVIParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
gBVIParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
if (force.getNonbondedMethod() != GBVISoftcoreForce::NoCutoff){
gBVIParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
}
if (force.getNonbondedMethod() == GBVISoftcoreForce::CutoffPeriodic) {
Vec3 boxVectors[3];
system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
RealOpenMM periodicBoxSize[3];
periodicBoxSize[0] = (RealOpenMM) boxVectors[0][0];
periodicBoxSize[1] = (RealOpenMM) boxVectors[1][1];
periodicBoxSize[2] = (RealOpenMM) boxVectors[2][2];
gBVIParameters->setPeriodic(periodicBoxSize);
}
gbviSoftcore = new CpuGBVISoftcore(gBVIParameters);
}
double ReferenceFreeEnergyCalcGBVISoftcoreForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealOpenMM> bornRadii(context.getSystem().getNumParticles());
gbviSoftcore->computeBornRadii(posData, bornRadii );
if (includeForces) {
vector<RealVec>& forceData = extractForces(context);
gbviSoftcore->computeBornForces(bornRadii, posData, &charges[0], forceData);
}
RealOpenMM energy = 0.0;
if (includeEnergy)
energy = gbviSoftcore->computeBornEnergy(bornRadii, posData, &charges[0]);
return static_cast<double>(energy);
}
plugins/freeEnergy/platforms/reference/src/ReferenceFreeEnergyKernels.h deleted 100644 → 0
View file @ c775bd19
#ifndef OPENMM_REFERENCE_FREE_ENERGY_KERNELS_H_
#define OPENMM_REFERENCE_FREE_ENERGY_KERNELS_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ReferencePlatform.h"
#include "openmm/freeEnergyKernels.h"
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include "SimTKReference/ReferenceNeighborList.h"
#include "gbsa/CpuGBVISoftcore.h"
#include "gbsa/CpuObcSoftcore.h"
namespace OpenMM {
/**
* This kernel is invoked by NonbondedSoftcoreForce to calculate the forces acting on the system.
*/
class ReferenceFreeEnergyCalcNonbondedSoftcoreForceKernel : public CalcNonbondedSoftcoreForceKernel {
public:
ReferenceFreeEnergyCalcNonbondedSoftcoreForceKernel(std::string name, const Platform& platform) : CalcNonbondedSoftcoreForceKernel(name, platform) {
}
~ReferenceFreeEnergyCalcNonbondedSoftcoreForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the NonbondedSoftcoreForce this kernel will be used for
*/
void initialize(const System& system, const NonbondedSoftcoreForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
private:
int numParticles, num14;
int **exclusionArray, **bonded14IndexArray;
RealOpenMM **particleParamArray, **bonded14ParamArray;
RealOpenMM nonbondedCutoff, rfDielectric, ewaldAlpha;
RealVec periodicBoxSize;
int kmax[3];
std::vector<std::set<int> > exclusions;
NonbondedSoftcoreMethod nonbondedMethod;
NeighborList* neighborList;
};
/**
* This kernel is invoked by GBSAOBCForce to calculate the forces acting on the system.
*/
class ReferenceFreeEnergyCalcGBSAOBCSoftcoreForceKernel : public CalcGBSAOBCSoftcoreForceKernel {
public:
ReferenceFreeEnergyCalcGBSAOBCSoftcoreForceKernel(std::string name, const Platform& platform) : CalcGBSAOBCSoftcoreForceKernel(name, platform) {
}
~ReferenceFreeEnergyCalcGBSAOBCSoftcoreForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the GBSAOBCSoftcoreForce this kernel will be used for
*/
void initialize(const System& system, const GBSAOBCSoftcoreForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
private:
CpuObcSoftcore* obc;
std::vector<RealOpenMM> charges;
};
/**
* This kernel is invoked by GBVISoftcoreForce to calculate the forces acting on the system.
*/
class ReferenceFreeEnergyCalcGBVISoftcoreForceKernel : public CalcGBVISoftcoreForceKernel {
public:
ReferenceFreeEnergyCalcGBVISoftcoreForceKernel(std::string name, const Platform& platform) : CalcGBVISoftcoreForceKernel(name, platform) {
}
~ReferenceFreeEnergyCalcGBVISoftcoreForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the GBVISoftcoreForce this kernel will be used for
* @param scaled radii the scaled radii (Eq. 5 of Labute paper)
*/
void initialize(const System& system, const GBVISoftcoreForce& force, const std::vector<double> & scaledRadii);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
private:
CpuGBVISoftcore* gbviSoftcore;
std::vector<RealOpenMM> charges;
};
} // namespace OpenMM
#endif /*OPENMM_REFERENCE_FREE_ENERGY_KERNELS_H_*/
plugins/freeEnergy/platforms/reference/src/SimTKReference/ReferenceFreeEnergyLJCoulomb14Softcore.cpp deleted 100644 → 0
View file @ c775bd19
/* 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 "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceFreeEnergyLJCoulomb14Softcore.h"
#include "ReferenceForce.h"
using std::vector;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
ReferenceFreeEnergyLJCoulomb14Softcore constructor
--------------------------------------------------------------------------------------- */
ReferenceFreeEnergyLJCoulomb14Softcore::ReferenceFreeEnergyLJCoulomb14Softcore( ) : cutoff(false) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceFreeEnergyLJCoulomb14Softcore::ReferenceFreeEnergyLJCoulomb14Softcore";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
ReferenceFreeEnergyLJCoulomb14Softcore destructor
--------------------------------------------------------------------------------------- */
ReferenceFreeEnergyLJCoulomb14Softcore::~ReferenceFreeEnergyLJCoulomb14Softcore( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceFreeEnergyLJCoulomb14Softcore::~ReferenceFreeEnergyLJCoulomb14Softcore";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
@param distance the cutoff distance
@param solventDielectric the dielectric constant of the bulk solvent
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulomb14Softcore::setUseCutoff( RealOpenMM distance, RealOpenMM solventDielectric ) {
cutoff = true;
cutoffDistance = distance;
krf = pow(cutoffDistance, (RealOpenMM)-3.0)*(solventDielectric-1.0f)/(2.0f*solventDielectric+1.0f);
crf = (1.0f/cutoffDistance)*(3.0f*solventDielectric)/(2.0f*solventDielectric+1.0f);
}
/**---------------------------------------------------------------------------------------
Calculate parameters for LJ 1-4 ixn
@param c6 c6
@param c12 c12
@param q1 q1 charge atom 1
@param q2 q2 charge atom 2
@param epsfac epsfac ????????????
@param parameters output parameters:
parameter[0]= c6*c6/c12
parameter[1]= (c12/c6)**1/6
parameter[2]= epsfactor*q1*q2
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulomb14Softcore::getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM q2, RealOpenMM epsfac,
RealOpenMM* parameters ) const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM six = 6.0;
static const RealOpenMM oneSixth = one/six;
// ---------------------------------------------------------------------------------------
if( c12 <= zero ){
parameters[0] = one;
parameters[1] = zero;
} else {
parameters[0] = (c6*c6)/c12;
parameters[1] = POW( (c12/c6), oneSixth );
}
parameters[2] = epsfac*q1*q2;
}
/**---------------------------------------------------------------------------------------
Calculate LJ 1-4 ixn
@param atomIndices atom indices of 4 atoms in bond
@param atomCoordinates atom coordinates
@param parameters three parameters:
parameters[0]= (c12/c6)**1/6 (sigma)
parameters[1]= c6*c6/c12 (4*epsilon)
parameters[2]= epsfac*q1*q2
@param forces force array (forces added to current values)
@param totalEnergy if not null, the energy will be added to this
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulomb14Softcore::calculateBondIxn( int* atomIndices, vector<RealVec>& atomCoordinates,
RealOpenMM* parameters, vector<RealVec>& forces,
RealOpenMM* totalEnergy ) const {
static const std::string methodName = "\nReferenceFreeEnergyLJCoulomb14Softcore::calculateBondIxn";
// 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;
// number of parameters
static const int numberOfParameters = 3;
// debug flag
static const int debug = 0;
static const int LastAtomIndex = 2;
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
// ---------------------------------------------------------------------------------------
// get deltaR, R2, and R between 2 atoms
int atomAIndex = atomIndices[0];
int atomBIndex = atomIndices[1];
ReferenceForce::getDeltaR( atomCoordinates[atomBIndex], atomCoordinates[atomAIndex], deltaR[0] );
if (cutoff && deltaR[0][ReferenceForce::RIndex] > cutoffDistance)
return;
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM sig = parameters[0];
RealOpenMM eps = parameters[1];
RealOpenMM minSoftCoreLJLambda = parameters[3];
RealOpenMM energy = zero;
RealOpenMM dEdR = zero;
if( minSoftCoreLJLambda < one ){
calculateOneSoftCoreLJ14Ixn( deltaR[0][ReferenceForce::RIndex], sig, eps, minSoftCoreLJLambda, &dEdR, &energy );
} else {
calculateOneLJ14Ixn( inverseR, sig, eps, &dEdR, &energy );
}
if (cutoff)
dEdR += parameters[2]*(inverseR-2.0f*krf*r2);
else
dEdR += parameters[2]*inverseR;
dEdR *= inverseR*inverseR;
// accumulate forces
for( int ii = 0; ii < 3; ii++ ){
RealOpenMM force = dEdR*deltaR[0][ii];
forces[atomAIndex][ii] += force;
forces[atomBIndex][ii] -= force;
}
if (cutoff)
energy += parameters[2]*(inverseR+krf*r2-crf);
else
energy += parameters[2]*inverseR;
// accumulate energies
if (totalEnergy != NULL)
*totalEnergy += energy;
}
/**---------------------------------------------------------------------------------------
Calculate LJ pair ixn between two atoms
@param inverseR 1/r
@param sig sigma
@param eps epsilon
@param dEdR output force factor
@param energy LJ energy
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulomb14Softcore::calculateOneLJ14Ixn( RealOpenMM inverseR, RealOpenMM sig, RealOpenMM eps,
RealOpenMM* dEdR, RealOpenMM* energy ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "\nReferenceLJ14CoulombIxn::calculateOneLJIxn";
// ---------------------------------------------------------------------------------------
// 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 six = 6.0;
static const RealOpenMM twelve = 12.0;
RealOpenMM sig2 = inverseR*sig;
sig2 *= sig2;
RealOpenMM sig6 = sig2*sig2*sig2;
*dEdR = eps*( twelve*sig6 - six )*sig6;
*energy += eps*(sig6-one)*sig6;
}
/**---------------------------------------------------------------------------------------
Calculate softcore LJ pair ixn between two atoms
@param r r
@param sig sigma
@param eps epsilon
@param lambda lambda
@param dEdR output force factor
@param energy LJ energy
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulomb14Softcore::calculateOneSoftCoreLJ14Ixn( RealOpenMM r, RealOpenMM sig, RealOpenMM eps,
RealOpenMM lambda,
RealOpenMM* dEdR, RealOpenMM* energy ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "\nReferenceFreeEnergyLJCoulomb14Softcore::calculateOneSoftCoreLJ14Ixn";
// ---------------------------------------------------------------------------------------
// 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 six = 6.0;
static const RealOpenMM twelve = 12.0;
static const RealOpenMM alphaLJ = 0.5;
#if 0
RealOpenMM dEdROrig = 0.0;
RealOpenMM E_Orig = 0.0;
static int maxPrint = 0;
calculateOneLJIxn( one/r, sig, eps, &dEdROrig, &E_Orig );
#endif
// soft-core LJ energy = lambda*4*eps*[ 1/{alphaLJ*(1-lambda) + (r/sig)**6}**2 - 1/{alphaLJ*(1-lambda) + (r/sig)**6} ]
eps *= lambda;
RealOpenMM sig2 = r/sig;
sig2 *= sig2;
RealOpenMM sig6 = sig2*sig2*sig2;
RealOpenMM softcoreLJTerm = alphaLJ*(one - lambda) + sig6;
RealOpenMM softcoreLJInv = one/softcoreLJTerm;
RealOpenMM softcoreLJInv2 = softcoreLJInv*softcoreLJInv;
*dEdR = eps*softcoreLJInv2*( twelve*softcoreLJInv - six )*sig6;
*energy += eps*softcoreLJInv*( softcoreLJInv - one );
#if 0
if( maxPrint++ < 5 ){
printf( "r=%14.6e sig=%14.6e eps=%14.6e lambda=%14.6e de[%14.6e %14.6e] e[%14.6e %14.6e] %14.6e %14.6e\n",
r, sig, eps/lambda, lambda, dEdROrig, *dEdR, E_Orig, *energy, softcoreLJInv, sig6 );
}
#endif
}
plugins/freeEnergy/platforms/reference/src/SimTKReference/ReferenceFreeEnergyLJCoulomb14Softcore.h deleted 100644 → 0
View file @ c775bd19
/* 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.
*/
#ifndef __ReferenceFreeEnergyLJCoulomb14Softcore_H__
#define __ReferenceFreeEnergyLJCoulomb14Softcore_H__
#include "SimTKReference/ReferenceBondIxn.h"
// ---------------------------------------------------------------------------------------
class ReferenceFreeEnergyLJCoulomb14Softcore : public ReferenceBondIxn {
private:
bool cutoff;
RealOpenMM cutoffDistance;
RealOpenMM krf, crf;
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
ReferenceFreeEnergyLJCoulomb14Softcore( );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceFreeEnergyLJCoulomb14Softcore( );
/**---------------------------------------------------------------------------------------
Set the force to use a cutoff.
@param distance the cutoff distance
@param solventDielectric the dielectric constant of the bulk solvent
--------------------------------------------------------------------------------------- */
void setUseCutoff( RealOpenMM distance, RealOpenMM solventDielectric );
/**---------------------------------------------------------------------------------------
Calculate parameters for LJ 1-4 ixn
@param c6 c6
@param c12 c12
@param q1 q1 charge atom 1
@param q2 q2 charge atom 2
@param epsfac epsfac ????????????
@param parameters output parameters:
parameter[0]= c6*c6/c12
parameter[1]= (c12/c6)**1/6
parameter[2]= epsfactor*q1*q2
--------------------------------------------------------------------------------------- */
void getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM q2, RealOpenMM epsfac,
RealOpenMM* parameters ) const;
/**---------------------------------------------------------------------------------------
Calculate Ryckaert-Bellemans bond ixn
@param atomIndices atom indices of 4 atoms in bond
@param atomCoordinates atom coordinates
@param parameters six RB parameters
@param forces force array (forces added to current values)
@param totalEnergy if not null, the energy will be added to this
--------------------------------------------------------------------------------------- */
void calculateBondIxn( int* atomIndices, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM* parameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* totalEnergy ) const;
/**---------------------------------------------------------------------------------------
Calculate LJ pair ixn between two atoms
@param inverseR 1/r
@param sig sigma
@param eps epsilon
@param dEdR output force factor
@param energy LJ energy
--------------------------------------------------------------------------------------- */
void calculateOneLJ14Ixn( RealOpenMM inverseR, RealOpenMM sig, RealOpenMM eps,
RealOpenMM* dEdR, RealOpenMM* energy ) const;
/**---------------------------------------------------------------------------------------
Calculate softcore LJ pair ixn between two atoms
@param r r
@param sig sigma
@param eps epsilon
@param lambda lambda
@param dEdR output force factor
@param energy LJ energy
--------------------------------------------------------------------------------------- */
void calculateOneSoftCoreLJ14Ixn( RealOpenMM r, RealOpenMM sig, RealOpenMM eps,
RealOpenMM lambda, RealOpenMM* dEdR, RealOpenMM* energy ) const;
};
// ---------------------------------------------------------------------------------------
#endif // __ReferenceFreeEnergyLJCoulomb14Softcore_H__
plugins/freeEnergy/platforms/reference/src/SimTKReference/ReferenceFreeEnergyLJCoulombSoftcoreIxn.cpp deleted 100644 → 0
View file @ c775bd19
/* 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 "ReferenceFreeEnergyLJCoulombSoftcoreIxn.h"
#include "ReferenceForce.h"
using std::vector;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
ReferenceFreeEnergyLJCoulombSoftcoreIxn constructor
--------------------------------------------------------------------------------------- */
ReferenceFreeEnergyLJCoulombSoftcoreIxn::ReferenceFreeEnergyLJCoulombSoftcoreIxn( ) : cutoff(false), periodic(false), ewald(false), pme(false) {
}
/**---------------------------------------------------------------------------------------
ReferenceFreeEnergyLJCoulombSoftcoreIxn destructor
--------------------------------------------------------------------------------------- */
ReferenceFreeEnergyLJCoulombSoftcoreIxn::~ReferenceFreeEnergyLJCoulombSoftcoreIxn( ){
}
/**---------------------------------------------------------------------------------------
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
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric ) {
cutoff = true;
cutoffDistance = distance;
neighborList = &neighbors;
krf = pow(cutoffDistance, (RealOpenMM)-3.0)*(solventDielectric-1.0f)/(2.0f*solventDielectric+1.0f);
crf = (1.0f/cutoffDistance)*(3.0f*solventDielectric)/(2.0f*solventDielectric+1.0f);
}
/**---------------------------------------------------------------------------------------
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
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::setPeriodic( RealVec& 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];
}
/**---------------------------------------------------------------------------------------
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
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM epsfacSqrt,
RealOpenMM* parameters ) const {
// ---------------------------------------------------------------------------------------
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;
}
/**---------------------------------------------------------------------------------------
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 totalEnergy total energy
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePairIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
std::vector<RealOpenMM> de( numberOfAtoms, 0.0 );
if (cutoff) {
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, totalEnergy, de);
}
} else {
// allocate and initialize exclusion array
std::vector<int> exclusionIndices( numberOfAtoms, -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, totalEnergy, de);
}
}
}
}
}
/**---------------------------------------------------------------------------------------
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 totalEnergy total energy
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj,
vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, vector<RealVec>& forces,
RealOpenMM* totalEnergy, std::vector<RealOpenMM>& de ) const {
// ---------------------------------------------------------------------------------------
// 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;
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 eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
RealOpenMM minSoftCoreLJLambda = atomParameters[ii][SoftCoreLJLambdaIndex] < atomParameters[jj][SoftCoreLJLambdaIndex] ?
atomParameters[ii][SoftCoreLJLambdaIndex] : atomParameters[jj][SoftCoreLJLambdaIndex];
// LJ: use soft core LJ if lambda < 1
RealOpenMM energy = zero;
RealOpenMM dEdR;
RealOpenMM dEdRx;
RealOpenMM dEdCol;
if( minSoftCoreLJLambda < one ){
calculateOneSoftCoreLJIxn( deltaR[0][ReferenceForce::RIndex], sig, eps, minSoftCoreLJLambda, &dEdR, &energy );
} else {
calculateOneLJIxn( inverseR, sig, eps, &dEdR, &energy );
}
dEdRx = dEdR;
// Coulomb
if (cutoff){
dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
dEdCol = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
} else {
dEdR += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
dEdCol = atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
}
dEdCol *= inverseR*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;
}
// accumulate energies
if( totalEnergy ) {
if (cutoff)
energy += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf);
else
energy += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
if( totalEnergy )
*totalEnergy += energy;
}
}
/**---------------------------------------------------------------------------------------
Calculate LJ pair ixn between two atoms
@param inverseR 1/r
@param sig sigma
@param eps epsilon
@param dEdR output force factor
@param energy LJ energy
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneLJIxn( RealOpenMM inverseR, RealOpenMM sig, RealOpenMM eps,
RealOpenMM* dEdR, RealOpenMM* energy ) const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM six = 6.0;
static const RealOpenMM twelve = 12.0;
RealOpenMM sig2 = inverseR*sig;
sig2 *= sig2;
RealOpenMM sig6 = sig2*sig2*sig2;
*dEdR = eps*( twelve*sig6 - six )*sig6;
*energy += eps*(sig6-one)*sig6;
}
/**---------------------------------------------------------------------------------------
Calculate softcore LJ pair ixn between two atoms
@param r r
@param sig sigma
@param eps epsilon
@param lambda lambda
@param dEdR output force factor
@param energy LJ energy
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneSoftCoreLJIxn( RealOpenMM r, RealOpenMM sig, RealOpenMM eps,
RealOpenMM lambda,
RealOpenMM* dEdR, RealOpenMM* energy ) const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM four = 4.0;
static const RealOpenMM six = 6.0;
static const RealOpenMM twelve = 12.0;
static const RealOpenMM alphaLJ = 0.5;
// soft-core LJ energy = lambda*4*eps*[ 1/{alphaLJ*(1-lambda) + (r/sig)**6}**2 - 1/{alphaLJ*(1-lambda) + (r/sig)**6} ]
eps *= lambda;
RealOpenMM sig2 = r/sig;
sig2 *= sig2;
RealOpenMM sig6 = sig2*sig2*sig2;
RealOpenMM softcoreLJTerm = alphaLJ*(one - lambda) + sig6;
RealOpenMM softcoreLJInv = one/softcoreLJTerm;
RealOpenMM softcoreLJInv2 = softcoreLJInv*softcoreLJInv;
*dEdR = eps*softcoreLJInv2*( twelve*softcoreLJInv - six )*sig6;
*energy += eps*softcoreLJInv*( softcoreLJInv - one );
}
plugins/freeEnergy/platforms/reference/src/SimTKReference/ReferenceFreeEnergyLJCoulombSoftcoreIxn.h deleted 100644 → 0
View file @ c775bd19
#ifndef __ReferenceFreeEnergyLJCoulombSoftcoreIxn_H__
#define __ReferenceFreeEnergyLJCoulombSoftcoreIxn_H__
/* 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 "SimTKReference/ReferencePairIxn.h"
#include "SimTKReference/ReferenceNeighborList.h"
#include <vector>
// ---------------------------------------------------------------------------------------
class ReferenceFreeEnergyLJCoulombSoftcoreIxn : public ReferencePairIxn {
private:
bool cutoff;
bool periodic;
bool ewald;
bool pme;
const OpenMM::NeighborList* neighborList;
RealOpenMM periodicBoxSize[3];
RealOpenMM cutoffDistance;
RealOpenMM krf, crf;
int numRx, numRy, numRz;
RealOpenMM alphaEwald;
// parameter indices
static const int SigIndex = 0;
static const int EpsIndex = 1;
static const int QIndex = 2;
static const int SoftCoreLJLambdaIndex = 3;
/**---------------------------------------------------------------------------------------
Calculate LJ Coulomb pair ixn between two atoms
@param atom1 the index of the first atom
@param atom2 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 totalEnergy total energy
--------------------------------------------------------------------------------------- */
void calculateOneIxn( int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* totalEnergy, std::vector<RealOpenMM>& de ) const;
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
ReferenceFreeEnergyLJCoulombSoftcoreIxn( );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~ReferenceFreeEnergyLJCoulombSoftcoreIxn( );
/**---------------------------------------------------------------------------------------
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
--------------------------------------------------------------------------------------- */
void setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric );
/**---------------------------------------------------------------------------------------
Set the force to use periodic boundary conditions. This requires that a cutoff has
already 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
--------------------------------------------------------------------------------------- */
void setPeriodic( OpenMM::RealVec& boxSize );
/**---------------------------------------------------------------------------------------
Calculate parameters for ixn
@param c6 c6
@param c12 c12
@param q1 q1 charge atom
@param epsfacSqrt epsfacSqrt (what is this?)
@param parameters output parameters:
parameter[SigIndex] = sqrt(c6*c6/c12)
parameter[EpsIndex] = 0.5*( (c12/c6)**1/6 )
parameter[QIndex] = epsfactorSqrt*q1
--------------------------------------------------------------------------------------- */
void getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM epsfacSqrt, RealOpenMM* parameters ) const;
/**---------------------------------------------------------------------------------------
Calculate LJ Coulomb pair ixn
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters (charges, c6, c12, ...) 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 totalEnergy total energy
--------------------------------------------------------------------------------------- */
void calculatePairIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
private:
/**---------------------------------------------------------------------------------------
Calculate LJ pair ixn between two atoms
@param inverseR 1/r
@param sig sigma
@param eps epsilon
@param dEdR output force factor
@param energy LJ energy
--------------------------------------------------------------------------------------- */
void calculateOneLJIxn( RealOpenMM inverseR, RealOpenMM sig, RealOpenMM eps,
RealOpenMM* dEdR, RealOpenMM* energy) const;
/**---------------------------------------------------------------------------------------
Calculate softcore LJ pair ixn between two atoms
@param r r
@param sig sigma
@param eps epsilon
@param lambda lambda
@param dEdR output force factor
@param energy LJ energy
--------------------------------------------------------------------------------------- */
void calculateOneSoftCoreLJIxn( RealOpenMM r, RealOpenMM sig, RealOpenMM eps,
RealOpenMM lambda, RealOpenMM* dEdR, RealOpenMM* energy ) const;
};
// ---------------------------------------------------------------------------------------
#endif // __ReferenceFreeEnergyLJCoulombSoftcoreIxn_H__
plugins/freeEnergy/platforms/reference/src/gbsa/CpuGBVISoftcore.cpp deleted 100644 → 0
View file @ c775bd19
/* Portions copyright (c) 2006-2009 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 <stdio.h>
#include <math.h>
#include <sstream>
#include <vector>
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKReference/ReferenceForce.h"
#include "CpuGBVISoftcore.h"
using std::vector;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
CpuGBVISoftcore constructor
gbviParameters gbviParameters object
--------------------------------------------------------------------------------------- */
CpuGBVISoftcore::CpuGBVISoftcore( GBVISoftcoreParameters* gbviParameters ) : _gbviParameters(gbviParameters) {
_switchDeriviative.resize(_gbviParameters->getNumberOfAtoms());
}
/**---------------------------------------------------------------------------------------
CpuGBVISoftcore destructor
--------------------------------------------------------------------------------------- */
CpuGBVISoftcore::~CpuGBVISoftcore( ){
}
/**---------------------------------------------------------------------------------------
Get GBVISoftcoreParameters reference
@return GBVISoftcoreParameters reference
--------------------------------------------------------------------------------------- */
GBVISoftcoreParameters* CpuGBVISoftcore::getGBVISoftcoreParameters( void ) const {
return _gbviParameters;
}
/**---------------------------------------------------------------------------------------
Set GBVISoftcoreParameters reference
@param GBVISoftcoreParameters reference
--------------------------------------------------------------------------------------- */
void CpuGBVISoftcore::setGBVISoftcoreParameters( GBVISoftcoreParameters* gbviParameters ){
_gbviParameters = gbviParameters;
}
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _obcSoftcoreParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
RealOpenMMVector& CpuGBVISoftcore::getSwitchDeriviative( void ){
return _switchDeriviative;
}
/**---------------------------------------------------------------------------------------
Compute quintic spline value and associated derviative
@param x value to compute spline at
@param rl lower cutoff value
@param ru upper cutoff value
@param outValue value of spline at x
@param outDerivative value of derivative of spline at x
--------------------------------------------------------------------------------------- */
void CpuGBVISoftcore::quinticSpline( RealOpenMM x, RealOpenMM rl, RealOpenMM ru,
RealOpenMM* outValue, RealOpenMM* outDerivative ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM minusSix = static_cast<RealOpenMM>( -6.0 );
static const RealOpenMM minusTen = static_cast<RealOpenMM>( -10.0 );
static const RealOpenMM minusThirty = static_cast<RealOpenMM>( -30.0 );
static const RealOpenMM fifteen = static_cast<RealOpenMM>( 15.0 );
static const RealOpenMM sixty = static_cast<RealOpenMM>( 60.0 );
// ---------------------------------------------------------------------------------------
RealOpenMM numerator = x - rl;
RealOpenMM denominator = ru - rl;
RealOpenMM ratio = numerator/denominator;
RealOpenMM ratio2 = ratio*ratio;
RealOpenMM ratio3 = ratio2*ratio;
*outValue = one + ratio3*(minusTen + fifteen*ratio + minusSix*ratio2);
*outDerivative = ratio2*(minusThirty + sixty*ratio + minusThirty*ratio2)/denominator;
}
/**---------------------------------------------------------------------------------------
Compute Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
and quintic splice switching function
@param atomicRadius3 atomic radius cubed
@param bornSum Born sum (volume integral)
@param gbviParameters Gbvi parameters (parameters used in spline
QuinticLowerLimitFactor & QuinticUpperBornRadiusLimit)
@param bornRadius output Born radius
@param switchDeriviative output switching function deriviative
--------------------------------------------------------------------------------------- */
void CpuGBVISoftcore::computeBornRadiiUsingQuinticSpline( RealOpenMM atomicRadius3, RealOpenMM bornSum,
GBVISoftcoreParameters* gbviParameters,
RealOpenMM& bornRadius, RealOpenMM* switchDeriviative ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM minusOne = static_cast<RealOpenMM>( -1.0 );
static const RealOpenMM minusThree = static_cast<RealOpenMM>( -3.0 );
static const RealOpenMM oneEighth = static_cast<RealOpenMM>( 0.125 );
static const RealOpenMM minusOneThird = static_cast<RealOpenMM>( (-1.0/3.0) );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
// ---------------------------------------------------------------------------------------
// R = [ S(V)*(A - V) ]**(-1/3)
// S(V) = 1 V < L
// S(V) = qSpline + U/(A-V) L < V < A
// S(V) = U/(A-V) U < V
// dR/dr = (-1/3)*[ S(V)*(A - V) ]**(-4/3)*[ d{ S(V)*(A-V) }/dr
// d{ S(V)*(A-V) }/dr = (dV/dr)*[ (A-V)*dS/dV - S(V) ]
// (A - V)*dS/dV - S(V) = 0 - 1 V < L
// (A - V)*dS/dV - S(V) = (A-V)*d(qSpline) + (A-V)*U/(A-V)**2 - qSpline - U/(A-V)
// = (A-V)*d(qSpline) - qSpline L < V < A**(-3)
// (A - V)*dS/dV - S(V) = (A-V)*U*/(A-V)**2 - U/(A-V) = 0 U < V
RealOpenMM splineL = gbviParameters->getQuinticLowerLimitFactor()*atomicRadius3;
RealOpenMM sum;
if( bornSum > splineL ){
if( bornSum < atomicRadius3 ){
RealOpenMM splineValue, splineDerivative;
quinticSpline( bornSum, splineL, atomicRadius3, &splineValue, &splineDerivative );
sum = (atomicRadius3 - bornSum)*splineValue + gbviParameters->getQuinticUpperSplineLimit();
*switchDeriviative = splineValue - (atomicRadius3 - bornSum)*splineDerivative;
} else {
sum = gbviParameters->getQuinticUpperSplineLimit();
*switchDeriviative = zero;
}
} else {
sum = atomicRadius3 - bornSum;
*switchDeriviative = one;
}
bornRadius = POW( sum, minusOneThird );
}
/**---------------------------------------------------------------------------------------
Get Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
@param chain not used here
--------------------------------------------------------------------------------------- */
void CpuGBVISoftcore::computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMMVector& bornRadii ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM minusThree = static_cast<RealOpenMM>( -3.0 );
static const RealOpenMM oneEighth = static_cast<RealOpenMM>( 0.125 );
static const RealOpenMM minusOneThird = static_cast<RealOpenMM>( (-1.0/3.0) );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
// ---------------------------------------------------------------------------------------
GBVISoftcoreParameters* gbviParameters = getGBVISoftcoreParameters();
int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMMVector& atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMMVector& scaledRadii = gbviParameters->getScaledRadii();
const RealOpenMMVector& bornRadiusScaleFactors = gbviParameters->getBornRadiusScaleFactors();
RealOpenMMVector& switchDeriviative = getSwitchDeriviative();
// ---------------------------------------------------------------------------------------
// calculate Born radii
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM sum = zero;
// sum over volumes
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (_gbviParameters->getUseCutoff() && r > _gbviParameters->getCutoffDistance())
continue;
sum += bornRadiusScaleFactors[atomJ]*CpuGBVISoftcore::getVolume( r, radiusI, scaledRadii[atomJ] );
}
}
RealOpenMM atomicRadius3 = POW( radiusI, minusThree );
if( _gbviParameters->getBornRadiusScalingSoftcoreMethod() == GBVISoftcoreParameters::NoScaling ){
sum = atomicRadius3 - sum;
bornRadii[atomI] = POW( sum, minusOneThird );
switchDeriviative[atomI] = one;
} else if( _gbviParameters->getBornRadiusScalingSoftcoreMethod() == GBVISoftcoreParameters::QuinticSpline ){
RealOpenMM switchDeriviativeValue;
computeBornRadiiUsingQuinticSpline( atomicRadius3, sum, gbviParameters,
bornRadii[atomI], &switchDeriviativeValue );
switchDeriviative[atomI] = switchDeriviativeValue;
}
}
}
/**---------------------------------------------------------------------------------------
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return volume
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVISoftcore::getVolume( RealOpenMM r, RealOpenMM R, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM minusThree = static_cast<RealOpenMM>( -3.0 );
RealOpenMM diff = (S - R);
if( FABS( diff ) < r ){
RealOpenMM lowerBound = (R > (r - S)) ? R : (r - S);
return (CpuGBVISoftcore::getL( r, (r + S), S ) -
CpuGBVISoftcore::getL( r, lowerBound, S ));
} else if( r <= diff ){
return CpuGBVISoftcore::getL( r, (r + S), S ) -
CpuGBVISoftcore::getL( r, (r - S), S ) +
POW( R, minusThree );
} else {
return zero;
}
}
/**---------------------------------------------------------------------------------------
Get L (used in analytical solution for volume integrals)
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVISoftcore::getL( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM threeHalves = static_cast<RealOpenMM>( 1.5 );
static const RealOpenMM third = static_cast<RealOpenMM>( (1.0/3.0) );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
static const RealOpenMM eighth = static_cast<RealOpenMM>( 0.125 );
// ---------------------------------------------------------------------------------------
RealOpenMM rInv = one/r;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM diff2 = (r + S)*(r - S);
return (threeHalves*xInv2)*( (fourth*rInv) - (xInv*third) + (diff2*xInv2*eighth*rInv) );
}
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt r
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVISoftcore::dL_dr( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM threeHalves = static_cast<RealOpenMM>( 1.5 );
static const RealOpenMM threeEights = static_cast<RealOpenMM>( 0.375 );
static const RealOpenMM third = static_cast<RealOpenMM>( (1.0/3.0) );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
static const RealOpenMM eighth = static_cast<RealOpenMM>( 0.125 );
// ---------------------------------------------------------------------------------------
RealOpenMM rInv = one/r;
RealOpenMM rInv2 = rInv*rInv;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv3 = xInv2*xInv;
RealOpenMM diff2 = (r + S)*(r - S);
return ( (-threeHalves*xInv2*rInv2)*( fourth + eighth*diff2*xInv2 ) + threeEights*xInv3*xInv );
}
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt x
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVISoftcore::dL_dx( RealOpenMM r, RealOpenMM x, RealOpenMM S ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM threeHalvesM = static_cast<RealOpenMM>( -1.5 );
static const RealOpenMM third = static_cast<RealOpenMM>( (1.0/3.0) );
// ---------------------------------------------------------------------------------------
RealOpenMM rInv = one/r;
RealOpenMM xInv = one/x;
RealOpenMM xInv2 = xInv*xInv;
RealOpenMM xInv3 = xInv2*xInv;
RealOpenMM diff = (r + S)*(r - S);
return (threeHalvesM*xInv3)*( (half*rInv) - xInv + (half*diff*xInv2*rInv) );
}
/**---------------------------------------------------------------------------------------
Sgb function
@param t r*r*G_i*G_j
@return Sgb (top of p. 1694 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVISoftcore::Sgb( RealOpenMM t ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
// ---------------------------------------------------------------------------------------
return ( (t != zero) ? one/SQRT( (one + (fourth*EXP( -t ))/t) ) : zero);
}
/**---------------------------------------------------------------------------------------
Get GB/VI energy
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM CpuGBVISoftcore::computeBornEnergy( const vector<RealOpenMM>& bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM two = static_cast<RealOpenMM>( 2.0 );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
static const RealOpenMM four = static_cast<RealOpenMM>( 4.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
static const RealOpenMM eighth = static_cast<RealOpenMM>( 0.125 );
// ---------------------------------------------------------------------------------------
const GBVISoftcoreParameters* gbviParameters = getGBVISoftcoreParameters();
const RealOpenMM preFactor = gbviParameters->getElectricConstant();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMMVector& atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMMVector& gammaParameters = gbviParameters->getGammaParameters();
// ---------------------------------------------------------------------------------------
// Eq.2 of Labute paper [JCC 29 p. 1693-1698 2008]
// to minimze roundoff error sum cavityEnergy separately since in general much
// smaller than other contributions
RealOpenMM energy = zero;
RealOpenMM cavityEnergy = zero;
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM partialChargeI = partialCharges[atomI];
// self-energy term
RealOpenMM atomIEnergy = half*partialChargeI/bornRadii[atomI];
// cavity term
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
cavityEnergy += gammaParameters[atomI]*ratio*ratio*ratio;
for( int atomJ = atomI + 1; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM t = fourth*r2/(bornRadii[atomI]*bornRadii[atomJ]);
atomIEnergy += partialCharges[atomJ]*Sgb( t )/deltaR[ReferenceForce::RIndex];
}
energy += two*partialChargeI*atomIEnergy;
}
energy *= preFactor;
energy -= cavityEnergy;
RealOpenMM conversion = static_cast<RealOpenMM>(gbviParameters->getTau());
return (conversion*energy);
}
/**---------------------------------------------------------------------------------------
Get GB/VI forces
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
--------------------------------------------------------------------------------------- */
void CpuGBVISoftcore::computeBornForces( const vector<RealOpenMM>& bornRadii, vector<RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, vector<RealVec>& inputForces ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM two = static_cast<RealOpenMM>( 2.0 );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
static const RealOpenMM four = static_cast<RealOpenMM>( 4.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM oneThird = static_cast<RealOpenMM>( (1.0/3.0) );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
static const RealOpenMM eighth = static_cast<RealOpenMM>( 0.125 );
// ---------------------------------------------------------------------------------------
const GBVISoftcoreParameters* gbviParameters = getGBVISoftcoreParameters();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMMVector& atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMMVector& gammaParameters = gbviParameters->getGammaParameters();
const RealOpenMM preFactor = two*gbviParameters->getElectricConstant();
// ---------------------------------------------------------------------------------------
// set energy/forces to zero
RealOpenMMVector bornForces( numberOfAtoms, 0.0 );
std::vector<RealVec> forces( numberOfAtoms, RealVec(0.0, 0.0, 0.0) );
// ---------------------------------------------------------------------------------------
// first main loop
#undef TARGET
#define TARGET -1926
RealOpenMMVector bornForcesT( numberOfAtoms, 0.0 );
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// partial of polar term wrt Born radius
// and (dGpol/dr)(dr/dx)
RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
RealOpenMM D_ij = r2/(four*alpha2_ij);
RealOpenMM expTerm = EXP( -D_ij );
RealOpenMM denominator2 = r2 + alpha2_ij*expTerm;
RealOpenMM denominator = SQRT( denominator2 );
RealOpenMM Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
RealOpenMM dGpol_dr = -Gpol*( one - fourth*expTerm )/denominator2;
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*( one + D_ij )/denominator2;
if( atomI != atomJ ){
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
if( atomJ == TARGET ){
bornForcesT[atomI] = dGpol_dalpha2_ij*bornRadii[atomI];
}
deltaX *= dGpol_dr;
deltaY *= dGpol_dr;
deltaZ *= dGpol_dr;
forces[atomI][0] += deltaX;
forces[atomI][1] += deltaY;
forces[atomI][2] += deltaZ;
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
}
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
if( atomI == TARGET ){
bornForcesT[atomJ] = dGpol_dalpha2_ij*bornRadii[atomJ];
}
}
}
// second main loop: (dGpol/dBornRadius)(dBornRadius/dr)(dr/dx)
// dGpol/dBornRadius) = bornForces[]
// dBornRadius/dr = (1/3)*(bR**4)*(dV/dr)
const RealOpenMMVector& scaledRadii = gbviParameters->getScaledRadii();
RealOpenMMVector& switchDeriviative = getSwitchDeriviative();
const RealOpenMMVector& bornRadiusScaleFactors = gbviParameters->getBornRadiusScaleFactors();
if( 0 ){
double bF = 0.0;
int count = 0;
fprintf( stderr, "PdeRef %d\n", TARGET );
RealOpenMM tau = static_cast<RealOpenMM>(gbviParameters->getTau());
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
bF += bornForcesT[atomI];
if( fabs( bornForcesT[atomI] ) > 0.0 ){
count++;
fprintf( stderr, "%6d %6d bFAdd=%15.7e bR=%15.7e\n", atomI, count, tau*bornForcesT[atomI], bornRadii[atomI] );
}
}
RealOpenMM ratio = (atomicRadii[TARGET]/bornRadii[TARGET]);
RealOpenMM bF1 = bF + (three*gammaParameters[TARGET]*ratio*ratio*ratio)/bornRadii[TARGET];
RealOpenMM b2 = bornRadii[TARGET]*bornRadii[TARGET];
bF1 *= switchDeriviative[TARGET]*oneThird*b2*b2;
fprintf( stderr, "sumbF Ref %6d %15.7e %15.7e %15.7e %15.7e\n", TARGET, bF, tau*bF1, bornRadii[TARGET], switchDeriviative[TARGET] );
}
// ---------------------------------------------------------------------------------------
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM R = atomicRadii[atomI];
// partial of cavity term wrt Born radius
RealOpenMM ratio = (atomicRadii[atomI]/bornRadii[atomI]);
bornForces[atomI] += (three*gammaParameters[atomI]*ratio*ratio*ratio)/bornRadii[atomI];
RealOpenMM b2 = bornRadii[atomI]*bornRadii[atomI];
bornForces[atomI] *= switchDeriviative[atomI]*oneThird*b2*b2;
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
RealOpenMM deltaX = atomCoordinates[atomJ][0] - atomCoordinates[atomI][0];
RealOpenMM deltaY = atomCoordinates[atomJ][1] - atomCoordinates[atomI][1];
RealOpenMM deltaZ = atomCoordinates[atomJ][2] - atomCoordinates[atomI][2];
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_gbviParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _gbviParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_gbviParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _gbviParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
deltaX = deltaR[ReferenceForce::XIndex];
deltaY = deltaR[ReferenceForce::YIndex];
deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM r = SQRT( r2 );
RealOpenMM S = scaledRadii[atomJ];
RealOpenMM diff = (S - R);
// find dRb/dr, where Rb is the Born radius
RealOpenMM de = CpuGBVISoftcore::dL_dr( r, r+S, S ) + CpuGBVISoftcore::dL_dx( r, r+S, S );
if( FABS( diff ) < r ){
if( R > (r - S) ){
de -= CpuGBVISoftcore::dL_dr( r, R, S );
} else {
de -= ( CpuGBVISoftcore::dL_dr( r, (r-S), S ) + CpuGBVISoftcore::dL_dx( r, (r-S), S ) );
}
} else if( r < (S - R) ){
de -= ( CpuGBVISoftcore::dL_dr( r, r-S, S ) + CpuGBVISoftcore::dL_dx( r, r-S, S ) );
}
// de = (dG/dRb)(dRb/dr)
de *= bornRadiusScaleFactors[atomJ]*bornForces[atomI]/r;
deltaX *= de;
deltaY *= de;
deltaZ *= de;
forces[atomI][0] += deltaX;
forces[atomI][1] += deltaY;
forces[atomI][2] += deltaZ;
forces[atomJ][0] -= deltaX;
forces[atomJ][1] -= deltaY;
forces[atomJ][2] -= deltaZ;
}
}
}
//printGbvi( atomCoordinates, partialCharges, bornRadii,bornForces, forces, "Reference: GBVI softcore post loop2", stderr );
// apply prefactor tau = (1/diel_solute - 1/diel_solvent)
RealOpenMM conversion = static_cast<RealOpenMM>(gbviParameters->getTau());
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
inputForces[atomI][0] += conversion*forces[atomI][0];
inputForces[atomI][1] += conversion*forces[atomI][1];
inputForces[atomI][2] += conversion*forces[atomI][2];
}
}
/**---------------------------------------------------------------------------------------
Print GB/VI parameters, radii, forces, ...
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param bornRadii Born radii (may be empty)
@param bornForces Born forces (may be empty)
@param forces forces (may be empty)
@param idString id string (who is calling)
@param log log file
--------------------------------------------------------------------------------------- */
void CpuGBVISoftcore::printGbvi( const std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges,
const RealOpenMMVector& bornRadii,
const RealOpenMMVector& bornForces,
const std::vector<OpenMM::RealVec>& forces,
const std::string& idString, FILE* log ){
// ---------------------------------------------------------------------------------------
const GBVISoftcoreParameters* gbviParameters = getGBVISoftcoreParameters();
const int numberOfAtoms = gbviParameters->getNumberOfAtoms();
const RealOpenMMVector& atomicRadii = gbviParameters->getAtomicRadii();
const RealOpenMMVector& gammaParameters = gbviParameters->getGammaParameters();
const RealOpenMMVector& bornRadiusScaleFactors = gbviParameters->getBornRadiusScaleFactors();
const RealOpenMM preFactor = 2.0*gbviParameters->getElectricConstant();
// ---------------------------------------------------------------------------------------
const RealOpenMMVector& scaledRadii = gbviParameters->getScaledRadii();
const RealOpenMMVector& switchDeriviative = getSwitchDeriviative();
RealOpenMM tau = static_cast<RealOpenMM>(gbviParameters->getTau());
(void) fprintf( log, "Reference Gbvi %s atoms=%d\n", idString.c_str(), numberOfAtoms );
(void) fprintf( log, " tau %15.7e\n", tau );
(void) fprintf( log, " scaleMethod %d (QuinticEnum=%d)\n",
gbviParameters->getBornRadiusScalingSoftcoreMethod(), GBVISoftcoreParameters::QuinticSpline );
(void) fprintf( log, " preFactor %15.7e)\n", preFactor );
for( unsigned int atomI = 0; atomI < static_cast<unsigned int>(numberOfAtoms); atomI++ ){
(void) fprintf( log, "%6d r=%15.7e rSc=%15.7e swd=%15.7e lmda=%4.2f tau*gam=%15.7e q=%15.7e", atomI,
atomicRadii[atomI],
scaledRadii[atomI],
switchDeriviative[atomI],
bornRadiusScaleFactors[atomI],
tau*gammaParameters[atomI],
partialCharges[atomI] );
if( bornRadii.size() > atomI ){
(void) fprintf( log, " bR=%15.7e", bornRadii[atomI] );
}
if( bornForces.size() > atomI ){
(void) fprintf( log, " tau*bF=%15.7e", tau*bornForces[atomI] );
}
(void) fprintf( log, "\n" );
}
return;
}
/**---------------------------------------------------------------------------------------
Use double precision
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return volume
--------------------------------------------------------------------------------------- */
double CpuGBVISoftcore::getVolumeD( double r, double R, double S ){
// ---------------------------------------------------------------------------------------
static const double zero = 0.0;
static const double minusThree = -3.0;
double diff = (S - R);
if( fabs( diff ) < r ){
double lowerBound = (R > (r - S)) ? R : (r - S);
return (CpuGBVISoftcore::getLD( r, (r + S), S ) -
CpuGBVISoftcore::getLD( r, lowerBound, S ));
} else if( r < diff ){
return CpuGBVISoftcore::getLD( r, (r + S), S ) -
CpuGBVISoftcore::getLD( r, (r - S), S ) +
pow( R, minusThree );
} else {
return zero;
}
}
/**---------------------------------------------------------------------------------------
Use double precision
Get L (used in analytical solution for volume integrals)
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
double CpuGBVISoftcore::getLD( double r, double x, double S ){
// ---------------------------------------------------------------------------------------
static const double one = 1.0;
static const double threeHalves = 1.5;
static const double third = 1.0/3.0;
static const double fourth = 0.25;
static const double eighth = 0.125;
// ---------------------------------------------------------------------------------------
double rInv = one/r;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double diff2 = (r + S)*(r - S);
return (threeHalves*xInv)*( (xInv*fourth*rInv) - (xInv2*third) + (diff2*xInv3*eighth*rInv) );
}
/**---------------------------------------------------------------------------------------
Use double precision
Get partial derivative of L wrt r
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
double CpuGBVISoftcore::dL_drD( double r, double x, double S ){
// ---------------------------------------------------------------------------------------
static const double one = 1.0;
static const double threeHalves = 1.5;
static const double threeEights = 0.375;
static const double third = 1.0/3.0;
static const double fourth = 0.25;
static const double eighth = 0.125;
// ---------------------------------------------------------------------------------------
double rInv = one/r;
double rInv2 = rInv*rInv;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double diff2 = (r + S)*(r - S);
return ( (-threeHalves*xInv2*rInv2)*( fourth + eighth*diff2*xInv2 ) + threeEights*xInv3*xInv );
}
/**---------------------------------------------------------------------------------------
Use double precision
Get partial derivative of L wrt x
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
double CpuGBVISoftcore::dL_dxD( double r, double x, double S ){
// ---------------------------------------------------------------------------------------
static const double one = 1.0;
static const double half = 0.5;
static const double threeHalvesM = -1.5;
static const double third = 1.0/3.0;
// ---------------------------------------------------------------------------------------
double rInv = one/r;
double xInv = one/x;
double xInv2 = xInv*xInv;
double xInv3 = xInv2*xInv;
double diff = (r + S)*(r - S);
return (threeHalvesM*xInv3)*( (half*rInv) - xInv + (half*diff*xInv2*rInv) );
}
plugins/freeEnergy/platforms/reference/src/gbsa/CpuGBVISoftcore.h deleted 100755 → 0
View file @ c775bd19
/* 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.
*/
#ifndef __CpuGBVISoftcore_H__
#define __CpuGBVISoftcore_H__
#include "GBVISoftcoreParameters.h"
// ---------------------------------------------------------------------------------------
class CpuGBVISoftcore {
private:
// GB/VI parameters
GBVISoftcoreParameters* _gbviParameters;
// vector containing switching function derivative
RealOpenMMVector _switchDeriviative;
public:
/**---------------------------------------------------------------------------------------
Constructor
@param implicitSolventParameters ImplicitSolventParameters reference
@return CpuImplicitSolvent object
--------------------------------------------------------------------------------------- */
CpuGBVISoftcore( GBVISoftcoreParameters* gbviParameters );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~CpuGBVISoftcore( );
/**---------------------------------------------------------------------------------------
Return GBVISoftcoreParameters
@return GBVISoftcoreParameters
--------------------------------------------------------------------------------------- */
GBVISoftcoreParameters* getGBVISoftcoreParameters( void ) const;
/**---------------------------------------------------------------------------------------
Set ImplicitSolventParameters
@param ImplicitSolventParameters
--------------------------------------------------------------------------------------- */
void setGBVISoftcoreParameters( GBVISoftcoreParameters* gbviParameters );
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _implicitSolventParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
RealOpenMMVector& getSwitchDeriviative( void );
/**---------------------------------------------------------------------------------------
Get Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
@param switchDeriviative derivative of switch function
--------------------------------------------------------------------------------------- */
void computeBornRadii( std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMMVector& bornRadii );
/**---------------------------------------------------------------------------------------
Get state
title title (optional)
@return state string
--------------------------------------------------------------------------------------- */
std::string getStateString( const char* title ) const;
/**---------------------------------------------------------------------------------------
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return volume
--------------------------------------------------------------------------------------- */
static RealOpenMM getVolume( RealOpenMM r, RealOpenMM R, RealOpenMM S );
/**---------------------------------------------------------------------------------------
Get L (analytical solution for volume integrals)
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static RealOpenMM getL( RealOpenMM r, RealOpenMM x, RealOpenMM S );
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt r
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static RealOpenMM dL_dr( RealOpenMM r, RealOpenMM x, RealOpenMM S );
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt x
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static RealOpenMM dL_dx( RealOpenMM r, RealOpenMM x, RealOpenMM S );
/**---------------------------------------------------------------------------------------
Sgb function
@param t r*r*G_i*G_j
@return Sgb (p. 1694 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static RealOpenMM Sgb( RealOpenMM t );
/**---------------------------------------------------------------------------------------
Get GB/VI energy
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM computeBornEnergy( const RealOpenMMVector& bornRadii, std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges );
/**---------------------------------------------------------------------------------------
Get GB/VI forces
@param bornRadii Born radii
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces output forces
--------------------------------------------------------------------------------------- */
void computeBornForces( const RealOpenMMVector& bornRadii, std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMM* partialCharges, std::vector<OpenMM::RealVec>& inputForces );
/**---------------------------------------------------------------------------------------
Get volume Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return volume
--------------------------------------------------------------------------------------- */
static double getVolumeD( double r, double R, double S );
/**---------------------------------------------------------------------------------------
Get L (analytical solution for volume integrals)
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return L value (Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static double getLD( double r, double x, double S );
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt r
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static double dL_drD( double r, double x, double S );
/**---------------------------------------------------------------------------------------
Get partial derivative of L wrt x
@param r distance between atoms i & j
@param R atomic radius
@param S scaled atomic radius
@return partial derivative based on Eq. 4 of Labute paper [JCC 29 p. 1693-1698 2008])
--------------------------------------------------------------------------------------- */
static double dL_dxD( double r, double x, double S );
/**---------------------------------------------------------------------------------------
Compute quintic spline value and associated derviative
@param x value to compute spline at
@param rl lower cutoff value
@param ru upper cutoff value
@param outValue value of spline at x
@param outDerivative value of derivative of spline at x
--------------------------------------------------------------------------------------- */
void quinticSpline( RealOpenMM x, RealOpenMM rl, RealOpenMM ru, RealOpenMM* outValue, RealOpenMM* outDerivative );
/**---------------------------------------------------------------------------------------
Compute Born radii based on Eq. 3 of Labute paper [JCC 29 p. 1693-1698 2008])
and quintic splice switching function
@param atomicRadius3 atomic radius cubed
@param bornSum Born sum (volume integral)
@param gbviParameters Gbvi parameters (parameters used in spline
QuinticLowerLimitFactor & QuinticUpperLimit)
@param bornRadius output Born radius
@param switchDeriviative output switching function deriviative
--------------------------------------------------------------------------------------- */
void computeBornRadiiUsingQuinticSpline( RealOpenMM atomicRadius3, RealOpenMM bornSum,
GBVISoftcoreParameters* gbviParameters,
RealOpenMM& bornRadius, RealOpenMM* switchDeriviative );
/**---------------------------------------------------------------------------------------
Print GB/VI parameters, radii, forces, ...
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param bornRadii Born radii (may be empty)
@param bornForces Born forces (may be empty)
@param forces forces (may be empty)
@param idString id string (who is calling)
@param log log file
--------------------------------------------------------------------------------------- */
void printGbvi( const std::vector<OpenMM::RealVec>& atomCoordinates, const RealOpenMM* partialCharges,
const RealOpenMMVector& bornRadii,
const RealOpenMMVector& bornForces,
const std::vector<OpenMM::RealVec>& forces,
const std::string& idString, FILE* log );
};
// ---------------------------------------------------------------------------------------
#endif // __CpuGBVISoftcore_H__
plugins/freeEnergy/platforms/reference/src/gbsa/CpuObcSoftcore.cpp deleted 100644 → 0
View file @ c775bd19
/* Portions copyright (c) 2006-2009 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 <stdlib.h>
#include <cmath>
#include <cstdio>
#include <vector>
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKReference/ReferenceForce.h"
#include "CpuObcSoftcore.h"
using std::vector;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
CpuObcSoftcore constructor
obcSoftcoreParameters obcSoftcoreParameters object
--------------------------------------------------------------------------------------- */
CpuObcSoftcore::CpuObcSoftcore( ObcSoftcoreParameters* obcSoftcoreParameters ) :
_obcSoftcoreParameters(obcSoftcoreParameters),
_includeAceApproximation(1) {
_obcChain.resize(_obcSoftcoreParameters->getNumberOfAtoms());
}
/**---------------------------------------------------------------------------------------
CpuObcSoftcore destructor
--------------------------------------------------------------------------------------- */
CpuObcSoftcore::~CpuObcSoftcore( ){
}
/**---------------------------------------------------------------------------------------
Get ObcSoftcoreParameters reference
@return ObcSoftcoreParameters reference
--------------------------------------------------------------------------------------- */
ObcSoftcoreParameters* CpuObcSoftcore::getObcSoftcoreParameters( void ) const {
return _obcSoftcoreParameters;
}
/**---------------------------------------------------------------------------------------
Set ObcSoftcoreParameters reference
@param ObcSoftcoreParameters reference
--------------------------------------------------------------------------------------- */
void CpuObcSoftcore::setObcSoftcoreParameters( ObcSoftcoreParameters* obcSoftcoreParameters ){
_obcSoftcoreParameters = obcSoftcoreParameters;
}
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _obcSoftcoreParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
RealOpenMMVector& CpuObcSoftcore::getObcChain( void ){
return _obcChain;
}
/**---------------------------------------------------------------------------------------
Return flag signalling whether AceApproximation for nonpolar term is to be included
@return flag
--------------------------------------------------------------------------------------- */
int CpuObcSoftcore::includeAceApproximation( void ) const {
return _includeAceApproximation;
}
/**---------------------------------------------------------------------------------------
Set flag indicating whether AceApproximation is to be included
@param includeAceApproximation new includeAceApproximation value
--------------------------------------------------------------------------------------- */
void CpuObcSoftcore::setIncludeAceApproximation( int includeAceApproximation ){
_includeAceApproximation = includeAceApproximation;
}
/**---------------------------------------------------------------------------------------
Calculation of Born radii based on papers:
J. Phys. Chem. 1996 100, 19824-19839 (HCT paper)
Proteins: Structure, Function, and Bioinformatcis 55:383-394 (2004) (OBC paper)
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
--------------------------------------------------------------------------------------- */
void CpuObcSoftcore::computeBornRadii( const vector<RealVec>& atomCoordinates, RealOpenMMVector& bornRadii ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM two = static_cast<RealOpenMM>( 2.0 );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
// ---------------------------------------------------------------------------------------
ObcSoftcoreParameters* obcSoftcoreParameters = getObcSoftcoreParameters();
int numberOfAtoms = obcSoftcoreParameters->getNumberOfAtoms();
const RealOpenMMVector& atomicRadii = obcSoftcoreParameters->getAtomicRadii();
const RealOpenMMVector& scaledRadiusFactor = obcSoftcoreParameters->getScaledRadiusFactors();
RealOpenMMVector& obcChain = getObcChain();
const RealOpenMMVector& nonPolarScaleFactors = obcSoftcoreParameters->getNonPolarScaleFactors();
RealOpenMM dielectricOffset = obcSoftcoreParameters->getDielectricOffset();
RealOpenMM alphaObc = obcSoftcoreParameters->getAlphaObc();
RealOpenMM betaObc = obcSoftcoreParameters->getBetaObc();
RealOpenMM gammaObc = obcSoftcoreParameters->getGammaObc();
// ---------------------------------------------------------------------------------------
// calculate Born radii
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
RealOpenMM radiusIInverse = one/offsetRadiusI;
RealOpenMM sum = zero;
// HCT code
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcSoftcoreParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcSoftcoreParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
RealOpenMM r = deltaR[ReferenceForce::RIndex];
if (_obcSoftcoreParameters->getUseCutoff() && r > _obcSoftcoreParameters->getCutoffDistance())
continue;
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
if( offsetRadiusI < rScaledRadiusJ ){
RealOpenMM rInverse = one/r;
RealOpenMM l_ij = offsetRadiusI > FABS( r - scaledRadiusJ ) ? offsetRadiusI : FABS( r - scaledRadiusJ );
l_ij = one/l_ij;
RealOpenMM u_ij = one/rScaledRadiusJ;
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
RealOpenMM ratio = LN( (u_ij/l_ij) );
RealOpenMM term = l_ij - u_ij + fourth*r*(u_ij2 - l_ij2) + ( half*rInverse*ratio) + (fourth*scaledRadiusJ*scaledRadiusJ*rInverse)*(l_ij2 - u_ij2);
// this case (atom i completely inside atom j) is not considered in the original paper
// Jay Ponder and the authors of Tinker recognized this and
// worked out the details
if( offsetRadiusI < (scaledRadiusJ - r) ){
term += two*( radiusIInverse - l_ij);
}
sum += nonPolarScaleFactors[atomJ]*term;
}
}
}
// OBC-specific code (Eqs. 6-8 in paper)
sum *= half*offsetRadiusI;
RealOpenMM sum2 = sum*sum;
RealOpenMM sum3 = sum*sum2;
RealOpenMM tanhSum = TANH( alphaObc*sum - betaObc*sum2 + gammaObc*sum3 );
bornRadii[atomI] = one/( one/offsetRadiusI - tanhSum/radiusI );
obcChain[atomI] = offsetRadiusI*( alphaObc - two*betaObc*sum + three*gammaObc*sum2 );
obcChain[atomI] = (one - tanhSum*tanhSum)*obcChain[atomI]/radiusI;
}
return;
}
/**---------------------------------------------------------------------------------------
Get nonpolar solvation force constribution via ACE approximation
@param obcSoftcoreParameters parameters
@param vdwRadii Vdw radii
@param bornRadii Born radii
@param energy energy (output): value is incremented from input value
@param forces forces: values are incremented from input values
--------------------------------------------------------------------------------------- */
void CpuObcSoftcore::computeAceNonPolarForce( const ObcSoftcoreParameters* obcSoftcoreParameters,
const vector<RealOpenMM>& bornRadii, RealOpenMM* energy,
vector<RealOpenMM>& forces ) const {
static const RealOpenMM zero = 0.0;
static const RealOpenMM six = 6.0;
// ---------------------------------------------------------------------------------------
// compute the nonpolar solvation via ACE approximation
const RealOpenMM probeRadius = obcSoftcoreParameters->getProbeRadius();
const RealOpenMM surfaceAreaFactor = obcSoftcoreParameters->getPi4Asolv();
const RealOpenMMVector& atomicRadii = obcSoftcoreParameters->getAtomicRadii();
const RealOpenMMVector& nonPolarScaleFactors = obcSoftcoreParameters->getNonPolarScaleFactors();
int numberOfAtoms = obcSoftcoreParameters->getNumberOfAtoms();
// the original ACE equation is based on Eq.2 of
// M. Schaefer, C. Bartels and M. Karplus, "Solution Conformations
// and Thermodynamics of Structured Peptides: Molecular Dynamics
// Simulation with an Implicit Solvation Model", J. Mol. Biol.,
// 284, 835-848 (1998) (ACE Method)
// The original equation includes the factor (atomicRadii[atomI]/bornRadii[atomI]) to the first power,
// whereas here the ratio is raised to the sixth power: (atomicRadii[atomI]/bornRadii[atomI])**6
// This modification was made by Jay Ponder and is based on observations that the change yields better correlations w/
// expected values. Jay did not think it was important enough to write up, so there is
// no paper to cite.
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
if( bornRadii[atomI] > zero ){
RealOpenMM r = atomicRadii[atomI] + probeRadius;
RealOpenMM ratio6 = POW( atomicRadii[atomI]/bornRadii[atomI], six );
RealOpenMM saTerm = nonPolarScaleFactors[atomI]*surfaceAreaFactor*r*r*ratio6;
*energy += saTerm;
forces[atomI] -= six*saTerm/bornRadii[atomI];
}
}
}
/**---------------------------------------------------------------------------------------
Get Obc Born energy and forces
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM CpuObcSoftcore::computeBornEnergyForces( vector<RealVec>& atomCoordinates,
const RealOpenMMVector& partialCharges,
vector<RealVec>& inputForces ){
// ---------------------------------------------------------------------------------------
static const RealOpenMM zero = static_cast<RealOpenMM>( 0.0 );
static const RealOpenMM one = static_cast<RealOpenMM>( 1.0 );
static const RealOpenMM two = static_cast<RealOpenMM>( 2.0 );
static const RealOpenMM three = static_cast<RealOpenMM>( 3.0 );
static const RealOpenMM four = static_cast<RealOpenMM>( 4.0 );
static const RealOpenMM half = static_cast<RealOpenMM>( 0.5 );
static const RealOpenMM fourth = static_cast<RealOpenMM>( 0.25 );
static const RealOpenMM eighth = static_cast<RealOpenMM>( 0.125 );
// ---------------------------------------------------------------------------------------
const ObcSoftcoreParameters* obcSoftcoreParameters = getObcSoftcoreParameters();
const int numberOfAtoms = obcSoftcoreParameters->getNumberOfAtoms();
// ---------------------------------------------------------------------------------------
// constants
const RealOpenMM preFactor = two*obcSoftcoreParameters->getElectricConstant()*(
(one/obcSoftcoreParameters->getSoluteDielectric()) -
(one/obcSoftcoreParameters->getSolventDielectric()) );
const RealOpenMM dielectricOffset = obcSoftcoreParameters->getDielectricOffset();
// ---------------------------------------------------------------------------------------
// compute Born radii
RealOpenMMVector bornRadii( numberOfAtoms );
computeBornRadii( atomCoordinates, bornRadii );
RealOpenMM obcEnergy = zero;
RealOpenMMVector bornForces( numberOfAtoms, 0.0 );
// ---------------------------------------------------------------------------------------
// N*( 8 + pow) ACE
// compute the nonpolar solvation via ACE approximation
if( includeAceApproximation() ){
computeAceNonPolarForce( obcSoftcoreParameters, bornRadii, &obcEnergy, bornForces );
}
const RealOpenMMVector& nonPolarScaleFactors = obcSoftcoreParameters->getNonPolarScaleFactors();
// ---------------------------------------------------------------------------------------
// first main loop
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
RealOpenMM partialChargeI = preFactor*partialCharges[atomI];
for( int atomJ = atomI; atomJ < numberOfAtoms; atomJ++ ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcSoftcoreParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcSoftcoreParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcSoftcoreParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcSoftcoreParameters->getCutoffDistance())
continue;
RealOpenMM r2 = deltaR[ReferenceForce::R2Index];
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM alpha2_ij = bornRadii[atomI]*bornRadii[atomJ];
RealOpenMM D_ij = r2/(four*alpha2_ij);
RealOpenMM expTerm = EXP( -D_ij );
RealOpenMM denominator2 = r2 + alpha2_ij*expTerm;
RealOpenMM denominator = SQRT( denominator2 );
// charges are assumed to be scaled on input so nonPolarScaleFactor is not needed
RealOpenMM Gpol = (partialChargeI*partialCharges[atomJ])/denominator;
RealOpenMM dGpol_dalpha2_ij = -half*Gpol*expTerm*( one + D_ij )/denominator2;
RealOpenMM dGpol_dr = -Gpol*( one - fourth*expTerm )/denominator2;
if( atomI != atomJ ){
bornForces[atomJ] += dGpol_dalpha2_ij*bornRadii[atomI];
deltaX *= dGpol_dr;
deltaY *= dGpol_dr;
deltaZ *= dGpol_dr;
inputForces[atomI][0] += deltaX;
inputForces[atomI][1] += deltaY;
inputForces[atomI][2] += deltaZ;
inputForces[atomJ][0] -= deltaX;
inputForces[atomJ][1] -= deltaY;
inputForces[atomJ][2] -= deltaZ;
} else {
Gpol *= half;
}
obcEnergy += Gpol;
bornForces[atomI] += dGpol_dalpha2_ij*bornRadii[atomJ];
}
}
// ---------------------------------------------------------------------------------------
// second main loop
RealOpenMMVector& obcChain = getObcChain();
const RealOpenMMVector& atomicRadii = obcSoftcoreParameters->getAtomicRadii();
const RealOpenMM alphaObc = obcSoftcoreParameters->getAlphaObc();
const RealOpenMM betaObc = obcSoftcoreParameters->getBetaObc();
const RealOpenMM gammaObc = obcSoftcoreParameters->getGammaObc();
const RealOpenMMVector& scaledRadiusFactor = obcSoftcoreParameters->getScaledRadiusFactors();
// compute factor that depends only on the outer loop index
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
bornForces[atomI] *= bornRadii[atomI]*bornRadii[atomI]*obcChain[atomI];
}
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
// radius w/ dielectric offset applied
RealOpenMM radiusI = atomicRadii[atomI];
RealOpenMM offsetRadiusI = radiusI - dielectricOffset;
for( int atomJ = 0; atomJ < numberOfAtoms; atomJ++ ){
if( atomJ != atomI ){
RealOpenMM deltaR[ReferenceForce::LastDeltaRIndex];
if (_obcSoftcoreParameters->getPeriodic())
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomI], atomCoordinates[atomJ], _obcSoftcoreParameters->getPeriodicBox(), deltaR );
else
ReferenceForce::getDeltaR( atomCoordinates[atomI], atomCoordinates[atomJ], deltaR );
if (_obcSoftcoreParameters->getUseCutoff() && deltaR[ReferenceForce::RIndex] > _obcSoftcoreParameters->getCutoffDistance())
continue;
RealOpenMM deltaX = deltaR[ReferenceForce::XIndex];
RealOpenMM deltaY = deltaR[ReferenceForce::YIndex];
RealOpenMM deltaZ = deltaR[ReferenceForce::ZIndex];
RealOpenMM r = deltaR[ReferenceForce::RIndex];
// radius w/ dielectric offset applied
RealOpenMM offsetRadiusJ = atomicRadii[atomJ] - dielectricOffset;
RealOpenMM scaledRadiusJ = offsetRadiusJ*scaledRadiusFactor[atomJ];
RealOpenMM scaledRadiusJ2 = scaledRadiusJ*scaledRadiusJ;
RealOpenMM rScaledRadiusJ = r + scaledRadiusJ;
// dL/dr & dU/dr are zero (this can be shown analytically)
// removed from calculation
if( offsetRadiusI < rScaledRadiusJ ){
RealOpenMM l_ij = offsetRadiusI > FABS( r - scaledRadiusJ ) ? offsetRadiusI : FABS( r - scaledRadiusJ );
l_ij = one/l_ij;
RealOpenMM u_ij = one/rScaledRadiusJ;
RealOpenMM l_ij2 = l_ij*l_ij;
RealOpenMM u_ij2 = u_ij*u_ij;
RealOpenMM rInverse = one/r;
RealOpenMM r2Inverse = rInverse*rInverse;
RealOpenMM t3 = eighth*(one + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + fourth*LN( u_ij/l_ij )*r2Inverse;
t3 *= nonPolarScaleFactors[atomJ];
RealOpenMM de = bornForces[atomI]*t3*rInverse;
deltaX *= de;
deltaY *= de;
deltaZ *= de;
inputForces[atomI][0] -= deltaX;
inputForces[atomI][1] -= deltaY;
inputForces[atomI][2] -= deltaZ;
inputForces[atomJ][0] += deltaX;
inputForces[atomJ][1] += deltaY;
inputForces[atomJ][2] += deltaZ;
}
}
}
}
return obcEnergy;;
}
plugins/freeEnergy/platforms/reference/src/gbsa/CpuObcSoftcore.h deleted 100755 → 0
View file @ c775bd19
/* 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.
*/
#ifndef __CpuObcSoftcore_H__
#define __CpuObcSoftcore_H__
#include "ObcSoftcoreParameters.h"
// ---------------------------------------------------------------------------------------
class CpuObcSoftcore {
private:
// GBSA/OBC parameters
ObcSoftcoreParameters* _obcSoftcoreParameters;
// arrays containing OBC chain derivative
RealOpenMMVector _obcChain;
// flag to signal whether ACE approximation
// is to be included
int _includeAceApproximation;
public:
/**---------------------------------------------------------------------------------------
Constructor
@param implicitSolventParameters ImplicitSolventParameters reference
@return CpuImplicitSolvent object
--------------------------------------------------------------------------------------- */
CpuObcSoftcore( ObcSoftcoreParameters* obcSoftcoreParameters );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~CpuObcSoftcore( );
/**---------------------------------------------------------------------------------------
Return ObcSoftcoreParameters
@return ObcSoftcoreParameters
--------------------------------------------------------------------------------------- */
ObcSoftcoreParameters* getObcSoftcoreParameters( void ) const;
/**---------------------------------------------------------------------------------------
Set ImplicitSolventParameters
@param ImplicitSolventParameters
--------------------------------------------------------------------------------------- */
void setObcSoftcoreParameters( ObcSoftcoreParameters* obcSoftcoreParameters );
/**---------------------------------------------------------------------------------------
Return flag signalling whether AceApproximation for nonpolar term is to be included
@return flag
--------------------------------------------------------------------------------------- */
int includeAceApproximation( void ) const;
/**---------------------------------------------------------------------------------------
Set flag indicating whether AceApproximation is to be included
@param includeAceApproximation new includeAceApproximation value
--------------------------------------------------------------------------------------- */
void setIncludeAceApproximation( int includeAceApproximation );
/**---------------------------------------------------------------------------------------
Return OBC chain derivative: size = _implicitSolventParameters->getNumberOfAtoms()
On first call, memory for array is allocated if not set
@return array
--------------------------------------------------------------------------------------- */
RealOpenMMVector& getObcChain( void );
/**---------------------------------------------------------------------------------------
Get Born radii based on OBC
@param atomCoordinates atomic coordinates
@param bornRadii output array of Born radii
@param obcChain output array of OBC chain derivative factors; if NULL,
then ignored
--------------------------------------------------------------------------------------- */
void computeBornRadii( const std::vector<OpenMM::RealVec>& atomCoordinates, RealOpenMMVector& bornRadii );
/**---------------------------------------------------------------------------------------
Get nonpolar solvation force constribution via ACE approximation
@param obcSoftcoreParameters parameters
@param vdwRadii Vdw radii
@param bornRadii Born radii
@param energy energy (output): value is incremented from input value
@param forces forces: values are incremented from input values
--------------------------------------------------------------------------------------- */
void computeAceNonPolarForce( const ObcSoftcoreParameters* obcSoftcoreParameters,
const RealOpenMMVector& bornRadii, RealOpenMM* energy,
RealOpenMMVector& forces ) const;
/**---------------------------------------------------------------------------------------
Get Born energy and forces based on OBC
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM computeBornEnergyForces( std::vector<OpenMM::RealVec>& atomCoordinates,
const RealOpenMMVector& partialCharges, std::vector<OpenMM::RealVec>& forces );
};
// ---------------------------------------------------------------------------------------
#endif // __CpuObcSoftcore_H__
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