Skip to content

GitLab

Commit f6352d88 authored by Michael Sherman's avatar Michael Sherman
Browse files

(1) Fixed more VC++ warnings; this is almost clean now. 

(2) Changed erfc.cpp to a header file which defines missing erf() and erfc() as static functions when using VC++ but otherwise does nothing.
(3) Replaced all nonstandard uses of real() and imag() methods as lvalues; instead work directly with the complex values.
(4) Fixed a bug in TestFindExclusions which was caught by VC++'s extensive STL debug mode self-examination.

There are still problems in hilbert.c and not all CUDA tests pass.
parent 351563b5
Hide whitespace changes
Inline Side-by-side
Showing with 104 additions and 96 deletions
platforms/cuda/include/CudaPlatform.h
View file @ f6352d88
......@@ -35,7 +35,7 @@
#include "openmm/Platform.h"
#include "CudaStreamFactory.h"
class _gpuContext;
struct _gpuContext;
namespace OpenMM {
......
platforms/cuda/src/CudaKernels.cpp
View file @ f6352d88
......@@ -77,7 +77,7 @@ static double calcEnergy(OpenMMContextImpl& context, System& system) {
double* posData = new double[positions.getSize()*3];
positions.saveToArray(posData);
vector<Vec3> pos(positions.getSize());
for (int i = 0; i < pos.size(); i++)
for (int i = 0; i < (int)pos.size(); i++)
pos[i] = Vec3(posData[3*i], posData[3*i+1], posData[3*i+2]);
delete[] posData;
refContext.setPositions(pos);
......@@ -164,7 +164,7 @@ void CudaCalcPeriodicTorsionForceKernel::initialize(const System& system, const
data.primaryKernel = this;
data.hasPeriodicTorsions = true;
numTorsions = force.getNumTorsions();
const float RadiansToDegrees = 180.0/3.14159265;
const float RadiansToDegrees = (float)(180.0/3.14159265);
vector<int> particle1(numTorsions);
vector<int> particle2(numTorsions);
vector<int> particle3(numTorsions);
......@@ -276,26 +276,26 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
c12[i] = (float) (4*depth*pow(radius, 12.0));
exclusionList[i].push_back(i);
}
for (int i = 0; i < exclusions.size(); i++) {
for (int i = 0; i < (int)exclusions.size(); i++) {
exclusionList[exclusions[i].first].push_back(exclusions[i].second);
exclusionList[exclusions[i].second].push_back(exclusions[i].first);
}
CudaNonbondedMethod method = NO_CUTOFF;
if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
gpuSetNonbondedCutoff(gpu, force.getCutoffDistance(), 78.3f);
gpuSetNonbondedCutoff(gpu, (float)force.getCutoffDistance(), 78.3f);
method = CUTOFF;
}
if (force.getNonbondedMethod() == NonbondedForce::CutoffPeriodic) {
Vec3 boxVectors[3];
force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
gpuSetPeriodicBoxSize(gpu, boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]);
gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
method = PERIODIC;
}
if (force.getNonbondedMethod() == NonbondedForce::Ewald) {
Vec3 boxVectors[3];
force.getPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
gpuSetPeriodicBoxSize(gpu, boxVectors[0][0], boxVectors[1][1], boxVectors[2][2]);
gpuSetPeriodicBoxSize(gpu, (float)boxVectors[0][0], (float)boxVectors[1][1], (float)boxVectors[2][2]);
method = EWALD;
}
......@@ -387,7 +387,7 @@ static void initializeIntegration(const System& system, CudaPlatform::PlatformDa
invMass1[i] = 1.0f/mass[particle1Index];
invMass2[i] = 1.0f/mass[particle2Index];
}
gpuSetConstraintParameters(gpu, particle1, particle2, distance, invMass1, invMass2, integrator.getConstraintTolerance(), 4);
gpuSetConstraintParameters(gpu, particle1, particle2, distance, invMass1, invMass2, (float)integrator.getConstraintTolerance(), 4);
// Initialize any terms that haven't already been handled by a Force.
......
platforms/cuda/src/kernels/gpu.cpp
View file @ f6352d88
......@@ -481,7 +481,7 @@ void gpuSetConstraintParameters(gpuContext gpu, const vector<int>& atom1, const
// Find how many constraints each atom is involved in.
vector<int> constraintCount(gpu->natoms, 0);
for (int i = 0; i < atom1.size(); i++) {
for (int i = 0; i < (int)atom1.size(); i++) {
constraintCount[atom1[i]]++;
constraintCount[atom2[i]]++;
}
......@@ -491,7 +491,7 @@ void gpuSetConstraintParameters(gpuContext gpu, const vector<int>& atom1, const
// connected to.
vector<map<int, float> > settleConstraints(gpu->natoms);
for (int i = 0; i < atom1.size(); i++) {
for (int i = 0; i < (int)atom1.size(); i++) {
if (constraintCount[atom1[i]] == 2 && constraintCount[atom2[i]] == 2) {
settleConstraints[atom1[i]][atom2[i]] = distance[i];
settleConstraints[atom2[i]][atom1[i]] = distance[i];
......@@ -501,7 +501,7 @@ void gpuSetConstraintParameters(gpuContext gpu, const vector<int>& atom1, const
// Now remove the ones that don't actually form closed loops of three atoms.
vector<int> settleClusters;
for (int i = 0; i < settleConstraints.size(); i++) {
for (int i = 0; i < (int)settleConstraints.size(); i++) {
if (settleConstraints[i].size() == 2) {
int partner1 = settleConstraints[i].begin()->first;
int partner2 = (++settleConstraints[i].begin())->first;
......@@ -524,7 +524,7 @@ void gpuSetConstraintParameters(gpuContext gpu, const vector<int>& atom1, const
gpu->psSettleParameter = psSettleParameter;
gpu->sim.pSettleParameter = psSettleParameter->_pDevStream[0];
gpu->sim.settleConstraints = settleClusters.size();
for (int i = 0; i < settleClusters.size(); i++) {
for (int i = 0; i < (int)settleClusters.size(); i++) {
int atom1 = settleClusters[i];
int atom2 = settleConstraints[atom1].begin()->first;
int atom3 = (++settleConstraints[atom1].begin())->first;
......@@ -570,7 +570,7 @@ void gpuSetConstraintParameters(gpuContext gpu, const vector<int>& atom1, const
map<int, ShakeCluster> clusters;
vector<bool> invalidForShake(gpu->natoms, false);
for (int i = 0; i < atom1.size(); i++) {
for (int i = 0; i < (int)atom1.size(); i++) {
if (isShakeAtom[atom1[i]])
continue; // This is being taken care of with SETTLE.
......@@ -669,7 +669,7 @@ void gpuSetConstraintParameters(gpuContext gpu, const vector<int>& atom1, const
// Find connected constraints for LINCS.
vector<int> lincsConstraints;
for (unsigned int i = 0; i < atom1.size(); i++)
for (unsigned i = 0; i < atom1.size(); i++)
if (!isShakeAtom[atom1[i]])
lincsConstraints.push_back(i);
int numLincs = (int) lincsConstraints.size();
......@@ -679,9 +679,9 @@ void gpuSetConstraintParameters(gpuContext gpu, const vector<int>& atom1, const
atomConstraints[atom2[lincsConstraints[i]]].push_back(i);
}
vector<vector<int> > linkedConstraints(numLincs);
for (unsigned int atom = 0; atom < atomConstraints.size(); atom++) {
for (unsigned int i = 0; i < atomConstraints[atom].size(); i++)
for (int j = 0; j < i; j++) {
for (unsigned atom = 0; atom < atomConstraints.size(); atom++) {
for (unsigned i = 0; i < atomConstraints[atom].size(); i++)
for (unsigned j = 0; j < i; j++) {
int c1 = atomConstraints[atom][i];
int c2 = atomConstraints[atom][j];
linkedConstraints[c1].push_back(c2);
......@@ -689,10 +689,10 @@ void gpuSetConstraintParameters(gpuContext gpu, const vector<int>& atom1, const
}
}
int maxLinks = 0;
for (unsigned int i = 0; i < linkedConstraints.size(); i++)
for (unsigned i = 0; i < linkedConstraints.size(); i++)
maxLinks = max(maxLinks, (int) linkedConstraints[i].size());
int maxAtomConstraints = 0;
for (unsigned int i = 0; i < atomConstraints.size(); i++)
for (unsigned i = 0; i < atomConstraints.size(); i++)
maxAtomConstraints = max(maxAtomConstraints, (int) atomConstraints[i].size());
// Fill in the CUDA streams.
......@@ -1491,10 +1491,10 @@ void gpuBuildExclusionList(gpuContext gpu)
// Mark which work units have exclusions.
for (int atom1 = 0; atom1 < gpu->exclusions.size(); ++atom1)
for (int atom1 = 0; atom1 < (int)gpu->exclusions.size(); ++atom1)
{
int x = atom1/grid;
for (int j = 0; j < gpu->exclusions[atom1].size(); ++j)
for (int j = 0; j < (int)gpu->exclusions[atom1].size(); ++j)
{
int atom2 = gpu->exclusions[atom1][j];
int y = atom2/grid;
......@@ -1502,10 +1502,10 @@ void gpuBuildExclusionList(gpuContext gpu)
pWorkList[cell] |= 1;
}
}
if (gpu->sim.paddedNumberOfAtoms > gpu->natoms)
if ((int)gpu->sim.paddedNumberOfAtoms > gpu->natoms)
{
int lastBlock = gpu->natoms/grid;
for (int i = 0; i < gpu->sim.workUnits; ++i)
for (int i = 0; i < (int)gpu->sim.workUnits; ++i)
{
int x = pWorkList[i]>>17;
int y = (pWorkList[i]>>2)&0x7FFF;
......@@ -1521,7 +1521,7 @@ void gpuBuildExclusionList(gpuContext gpu)
unsigned int* pExclusionIndex = psExclusionIndex->_pSysData;
gpu->sim.pExclusionIndex = psExclusionIndex->_pDevData;
int numWithExclusions = 0;
for (int i = 0; i < psExclusionIndex->_length; ++i)
for (int i = 0; i < (int)psExclusionIndex->_length; ++i)
if ((pWorkList[i]&1) == 1)
pExclusionIndex[i] = (numWithExclusions++)*grid;
......@@ -1531,13 +1531,13 @@ void gpuBuildExclusionList(gpuContext gpu)
gpu->psExclusion = psExclusion;
unsigned int* pExclusion = psExclusion->_pSysData;
gpu->sim.pExclusion = psExclusion->_pDevData;
for (int i = 0; i < psExclusion->_length; ++i)
for (int i = 0; i < (int)psExclusion->_length; ++i)
pExclusion[i] = 0xFFFFFFFF;
for (int atom1 = 0; atom1 < gpu->exclusions.size(); ++atom1)
for (int atom1 = 0; atom1 < (int)gpu->exclusions.size(); ++atom1)
{
int x = atom1/grid;
int offset1 = atom1-x*grid;
for (int j = 0; j < gpu->exclusions[atom1].size(); ++j)
for (int j = 0; j < (int)gpu->exclusions[atom1].size(); ++j)
{
int atom2 = gpu->exclusions[atom1][j];
int y = atom2/grid;
......@@ -1557,11 +1557,11 @@ void gpuBuildExclusionList(gpuContext gpu)
// Mark all interactions that involve a padding atom as being excluded.
for (int atom1 = gpu->natoms; atom1 < atoms; ++atom1)
for (int atom1 = gpu->natoms; atom1 < (int)atoms; ++atom1)
{
int x = atom1/grid;
int offset1 = atom1-x*grid;
for (int atom2 = 0; atom2 < atoms; ++atom2)
for (int atom2 = 0; atom2 < (int)atoms; ++atom2)
{
int y = atom2/grid;
int index = x*atoms+y*grid+offset1;
......@@ -1609,7 +1609,7 @@ static void tagAtomsInMolecule(int atom, int molecule, vector<int>& atomMolecule
// Recursively tag atoms as belonging to a particular molecule.
atomMolecule[atom] = molecule;
for (int i = 0; i < atomBonds[atom].size(); i++)
for (int i = 0; i < (int)atomBonds[atom].size(); i++)
if (atomMolecule[atomBonds[atom][i]] == -1)
tagAtomsInMolecule(atomBonds[atom][i], molecule, atomMolecule, atomBonds);
}
......@@ -1619,7 +1619,7 @@ static void findMoleculeGroups(gpuContext gpu)
// First make a list of constraints for future use.
vector<Constraint> constraints;
for (int i = 0; i < gpu->sim.ShakeConstraints; i++)
for (int i = 0; i < (int)gpu->sim.ShakeConstraints; i++)
{
int atom1 = (*gpu->psShakeID)[i].x;
int atom2 = (*gpu->psShakeID)[i].y;
......@@ -1632,7 +1632,7 @@ static void findMoleculeGroups(gpuContext gpu)
if (atom4 != -1)
constraints.push_back(Constraint(atom1, atom4, distance2));
}
for (int i = 0; i < gpu->sim.settleConstraints; i++)
for (int i = 0; i < (int)gpu->sim.settleConstraints; i++)
{
int atom1 = (*gpu->psSettleID)[i].x;
int atom2 = (*gpu->psSettleID)[i].y;
......@@ -1643,7 +1643,7 @@ static void findMoleculeGroups(gpuContext gpu)
constraints.push_back(Constraint(atom1, atom3, distance12*distance12));
constraints.push_back(Constraint(atom2, atom3, distance23*distance23));
}
for (int i = 0; i < gpu->sim.lincsConstraints; i++)
for (int i = 0; i < (int)gpu->sim.lincsConstraints; i++)
{
int atom1 = (*gpu->psLincsAtoms)[i].x;
int atom2 = (*gpu->psLincsAtoms)[i].y;
......@@ -1655,14 +1655,14 @@ static void findMoleculeGroups(gpuContext gpu)
int numAtoms = gpu->natoms;
vector<vector<int> > atomBonds(numAtoms);
for (int i = 0; i < gpu->sim.bonds; i++)
for (int i = 0; i < (int)gpu->sim.bonds; i++)
{
int atom1 = (*gpu->psBondID)[i].x;
int atom2 = (*gpu->psBondID)[i].y;
atomBonds[atom1].push_back(atom2);
atomBonds[atom2].push_back(atom1);
}
for (int i = 0; i < constraints.size(); i++)
for (int i = 0; i < (int)constraints.size(); i++)
{
int atom1 = constraints[i].atom1;
int atom2 = constraints[i].atom2;
......@@ -1686,27 +1686,27 @@ static void findMoleculeGroups(gpuContext gpu)
vector<Molecule> molecules(numMolecules);
for (int i = 0; i < numMolecules; i++)
molecules[i].atoms = atomIndices[i];
for (int i = 0; i < gpu->sim.bonds; i++)
for (int i = 0; i < (int)gpu->sim.bonds; i++)
{
int atom1 = (*gpu->psBondID)[i].x;
molecules[atomMolecule[atom1]].bonds.push_back(i);
}
for (int i = 0; i < gpu->sim.bond_angles; i++)
for (int i = 0; i < (int)gpu->sim.bond_angles; i++)
{
int atom1 = (*gpu->psBondAngleID1)[i].x;
molecules[atomMolecule[atom1]].angles.push_back(i);
}
for (int i = 0; i < gpu->sim.dihedrals; i++)
for (int i = 0; i < (int)gpu->sim.dihedrals; i++)
{
int atom1 = (*gpu->psDihedralID1)[i].x;
molecules[atomMolecule[atom1]].periodicTorsions.push_back(i);
}
for (int i = 0; i < gpu->sim.rb_dihedrals; i++)
for (int i = 0; i < (int)gpu->sim.rb_dihedrals; i++)
{
int atom1 = (*gpu->psRbDihedralID1)[i].x;
molecules[atomMolecule[atom1]].rbTorsions.push_back(i);
}
for (int i = 0; i < constraints.size(); i++)
for (int i = 0; i < (int)constraints.size(); i++)
{
molecules[atomMolecule[constraints[i].atom1]].constraints.push_back(i);
}
......@@ -1715,14 +1715,14 @@ static void findMoleculeGroups(gpuContext gpu)
vector<Molecule> uniqueMolecules;
vector<vector<int> > moleculeInstances;
for (int molIndex = 0; molIndex < molecules.size(); molIndex++)
for (int molIndex = 0; molIndex < (int)molecules.size(); molIndex++)
{
Molecule& mol = molecules[molIndex];
// See if it is identical to another molecule.
bool isNew = true;
for (int j = 0; j < uniqueMolecules.size() && isNew; j++)
for (int j = 0; j < (int)uniqueMolecules.size() && isNew; j++)
{
Molecule& mol2 = uniqueMolecules[j];
bool identical = true;
......@@ -1734,20 +1734,20 @@ static void findMoleculeGroups(gpuContext gpu)
float4* posq = gpu->psPosq4->_pSysData;
float4* velm = gpu->psVelm4->_pSysData;
float2* sigeps = gpu->psSigEps2->_pSysData;
for (int i = 0; i < mol.atoms.size() && identical; i++)
for (int i = 0; i < (int)mol.atoms.size() && identical; i++)
if (mol.atoms[i] != mol2.atoms[i]-atomOffset || posq[mol.atoms[i]].w != posq[mol2.atoms[i]].w ||
velm[mol.atoms[i]].w != velm[mol2.atoms[i]].w || sigeps[mol.atoms[i]].x != sigeps[mol2.atoms[i]].x ||
sigeps[mol.atoms[i]].y != sigeps[mol2.atoms[i]].y)
identical = false;
int4* bondID = gpu->psBondID->_pSysData;
float2* bondParam = gpu->psBondParameter->_pSysData;
for (int i = 0; i < mol.bonds.size() && identical; i++)
for (int i = 0; i < (int)mol.bonds.size() && identical; i++)
if (bondID[mol.bonds[i]].x != bondID[mol2.bonds[i]].x-atomOffset || bondID[mol.bonds[i]].y != bondID[mol2.bonds[i]].y-atomOffset ||
bondParam[mol.bonds[i]].x != bondParam[mol2.bonds[i]].x || bondParam[mol.bonds[i]].y != bondParam[mol2.bonds[i]].y)
identical = false;
int4* angleID = gpu->psBondAngleID1->_pSysData;
float2* angleParam = gpu->psBondAngleParameter->_pSysData;
for (int i = 0; i < mol.angles.size() && identical; i++)
for (int i = 0; i < (int)mol.angles.size() && identical; i++)
if (angleID[mol.angles[i]].x != angleID[mol2.angles[i]].x-atomOffset ||
angleID[mol.angles[i]].y != angleID[mol2.angles[i]].y-atomOffset ||
angleID[mol.angles[i]].z != angleID[mol2.angles[i]].z-atomOffset ||
......@@ -1756,7 +1756,7 @@ static void findMoleculeGroups(gpuContext gpu)
identical = false;
int4* periodicID = gpu->psDihedralID1->_pSysData;
float4* periodicParam = gpu->psDihedralParameter->_pSysData;
for (int i = 0; i < mol.periodicTorsions.size() && identical; i++)
for (int i = 0; i < (int)mol.periodicTorsions.size() && identical; i++)
if (periodicID[mol.periodicTorsions[i]].x != periodicID[mol2.periodicTorsions[i]].x-atomOffset ||
periodicID[mol.periodicTorsions[i]].y != periodicID[mol2.periodicTorsions[i]].y-atomOffset ||
periodicID[mol.periodicTorsions[i]].z != periodicID[mol2.periodicTorsions[i]].z-atomOffset ||
......@@ -1768,7 +1768,7 @@ static void findMoleculeGroups(gpuContext gpu)
int4* rbID = gpu->psRbDihedralID1->_pSysData;
float4* rbParam1 = gpu->psRbDihedralParameter1->_pSysData;
float2* rbParam2 = gpu->psRbDihedralParameter2->_pSysData;
for (int i = 0; i < mol.rbTorsions.size() && identical; i++)
for (int i = 0; i < (int)mol.rbTorsions.size() && identical; i++)
if (rbID[mol.rbTorsions[i]].x != rbID[mol2.rbTorsions[i]].x-atomOffset ||
rbID[mol.rbTorsions[i]].y != rbID[mol2.rbTorsions[i]].y-atomOffset ||
rbID[mol.rbTorsions[i]].z != rbID[mol2.rbTorsions[i]].z-atomOffset ||
......@@ -1780,7 +1780,7 @@ static void findMoleculeGroups(gpuContext gpu)
rbParam2[mol.rbTorsions[i]].x != rbParam2[mol2.rbTorsions[i]].x ||
rbParam2[mol.rbTorsions[i]].y != rbParam2[mol2.rbTorsions[i]].y)
identical = false;
for (int i = 0; i < mol.constraints.size() && identical; i++)
for (int i = 0; i < (int)mol.constraints.size() && identical; i++)
if (constraints[mol.constraints[i]].atom1 != constraints[mol2.constraints[i]].atom1-atomOffset ||
constraints[mol.constraints[i]].atom2 != constraints[mol2.constraints[i]].atom2-atomOffset ||
constraints[mol.constraints[i]].distance2 != constraints[mol2.constraints[i]].distance2)
......@@ -1799,12 +1799,12 @@ static void findMoleculeGroups(gpuContext gpu)
}
}
gpu->moleculeGroups.resize(moleculeInstances.size());
for (int i = 0; i < moleculeInstances.size(); i++)
for (int i = 0; i < (int)moleculeInstances.size(); i++)
{
gpu->moleculeGroups[i].instances = moleculeInstances[i];
vector<int>& atoms = uniqueMolecules[i].atoms;
gpu->moleculeGroups[i].atoms.resize(atoms.size());
for (int j = 0; j < atoms.size(); j++)
for (int j = 0; j < (int)atoms.size(); j++)
gpu->moleculeGroups[i].atoms[j] = atoms[j]-atoms[0];
}
}
......@@ -1852,7 +1852,7 @@ void gpuReorderAtoms(gpuContext gpu)
vector<int> originalIndex(numAtoms);
vector<float4> newPosq(numAtoms);
vector<float4> newVelm(numAtoms);
for (int group = 0; group < gpu->moleculeGroups.size(); group++)
for (int group = 0; group < (int)gpu->moleculeGroups.size(); group++)
{
// Find the center of each molecule.
......@@ -1865,7 +1865,7 @@ void gpuReorderAtoms(gpuContext gpu)
molPos[i].x = 0.0f;
molPos[i].y = 0.0f;
molPos[i].z = 0.0f;
for (int j = 0; j < atoms.size(); j++)
for (int j = 0; j < (int)atoms.size(); j++)
{
int atom = atoms[j]+mol.instances[i];
molPos[i].x += posq[atom].x;
......@@ -1890,7 +1890,7 @@ void gpuReorderAtoms(gpuContext gpu)
molPos[i].x -= dx;
molPos[i].y -= dy;
molPos[i].z -= dz;
for (int j = 0; j < atoms.size(); j++)
for (int j = 0; j < (int)atoms.size(); j++)
{
int atom = atoms[j]+mol.instances[i];
posq[atom].x -= dx;
......@@ -1906,9 +1906,9 @@ void gpuReorderAtoms(gpuContext gpu)
bool useHilbert = (numMolecules > 5000 || atoms.size() > 8); // For small systems, a simple zigzag curve works better than a Hilbert curve.
float binWidth;
if (useHilbert)
binWidth = max(max(maxx-minx, maxy-miny), maxz-minz)/255.0;
binWidth = (float)(max(max(maxx-minx, maxy-miny), maxz-minz)/255.0);
else
binWidth = 0.2*sqrt(gpu->sim.nonbondedCutoffSqr);
binWidth = (float)(0.2*sqrt(gpu->sim.nonbondedCutoffSqr));
int xbins = 1 + (int) ((maxx-minx)/binWidth);
int ybins = 1 + (int) ((maxy-miny)/binWidth);
vector<pair<int, int> > molBins(numMolecules);
......@@ -1942,7 +1942,7 @@ void gpuReorderAtoms(gpuContext gpu)
for (int i = 0; i < numMolecules; i++)
{
for (int j = 0; j < atoms.size(); j++)
for (int j = 0; j < (int)atoms.size(); j++)
{
int oldIndex = mol.instances[molBins[i].second]+atoms[j];
int newIndex = mol.instances[i]+atoms[j];
......
platforms/reference/src/SimTKReference/erfc.cpp → platforms/reference/src/SimTKReference/MSVC_erfc.h
View file @ f6352d88
#ifndef OPENMM_MSVC_ERFC_H_
#define OPENMM_MSVC_ERFC_H_
/*
* At least up to version 8 (VC++ 2005), Microsoft does not support the
* standard C99 erf() and erfc() functions. For now we're including these
* definitions for an MSVC compilation; if these are added later then
* the #ifdef below should change to compare _MSC_VER with a particular
* version level.
*/
#ifdef _MSC_VER
/***************************
* erf.cpp
......@@ -5,15 +18,15 @@
* written: 29-Jan-04
***************************/
#include <math.h>
#include <cmath>
static const double rel_error= 1E-12; //calculate 12 significant figures
//you can adjust rel_error to trade off between accuracy and speed
//but don't ask for > 15 figures (assuming usual 52 bit mantissa in a double)
static double erfc(double x);
double erf(double x)
static double erf(double x)
//erf(x) = 2/sqrt(pi)*integral(exp(-t^2),t,0,x)
// = 2/sqrt(pi)*[x - x^3/3 + x^5/5*2! - x^7/7*3! + ...]
// = 1-erfc(x)
......@@ -36,7 +49,7 @@ double erf(double x)
}
double erfc(double x)
static double erfc(double x)
//erfc(x) = 2/sqrt(pi)*integral(exp(-t^2),t,x,inf)
// = exp(-x^2)/sqrt(pi) * [1/x+ (1/2)/x+ (2/2)/x+ (3/2)/x+ (4/2)/x+ ...]
// = 1-erf(x)
......@@ -46,7 +59,8 @@ double erfc(double x)
if (fabs(x) < 2.2) {
return 1.0 - erf(x); //use series when fabs(x) < 2.2
}
if (signbit(x)) { //continued fraction only valid for x>0
// Don't look for x==0 here!
if (x < 0) { //continued fraction only valid for x>0
return 2.0 - erfc(-x);
}
double a=1, b=x; //last two convergent numerators
......@@ -66,3 +80,7 @@ double erfc(double x)
} while (fabs(q1-q2)/q2 > rel_error);
return one_sqrtpi*exp(-x*x)*q2;
}
#endif // _MSC_VER
#endif // OPENMM_MSVC_ERFC_H_
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
View file @ f6352d88
......@@ -32,7 +32,9 @@
#include "ReferenceLJCoulombIxn.h"
#include "ReferenceForce.h"
double erfc(double x);
// In case we're using some primitive version of Visual Studio this will
// make sure that erf() and erfc() are defined.
#include "MSVC_erfc.h"
/**---------------------------------------------------------------------------------------
......@@ -248,10 +250,10 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
RealOpenMM SQRT_PI = sqrt(PI);
RealOpenMM TWO_PI = 2.0 * PI;
static const RealOpenMM epsilon = 1.0;
static const RealOpenMM epsilon = 1.0;
static const RealOpenMM one = 1.0;
RealOpenMM recipCoeff = (RealOpenMM) 4.0*M_PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon;
RealOpenMM recipCoeff = (RealOpenMM)(4*M_PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
RealOpenMM selfEwaldEnergy = 0.0;
RealOpenMM realSpaceEwaldEnergy = 0.0;
......@@ -292,22 +294,16 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
}
for(i = 0; (i < numberOfAtoms); i++) {
for(m = 0; (m < 3); m++) {
EIR(0, i, m).real() = 1;
EIR(0, i, m).imag() = 0;
}
for(m = 0; (m < 3); m++)
EIR(0, i, m) = d_complex(1,0);
for(m=0; (m<3); m++) {
EIR(1, i, m).real() = cos(atomCoordinates[i][m]*recipBoxSize[m]);
EIR(1, i, m).imag() = sin(atomCoordinates[i][m]*recipBoxSize[m]);
}
for(m=0; (m<3); m++)
EIR(1, i, m) = d_complex(cos(atomCoordinates[i][m]*recipBoxSize[m]),
sin(atomCoordinates[i][m]*recipBoxSize[m]));
for(j=2; (j<kmax); j++) {
for(m=0; (m<3); m++) {
for(j=2; (j<kmax); j++)
for(m=0; (m<3); m++)
EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
}
}
}
// calculate reciprocal space energy and forces
......@@ -324,33 +320,25 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
RealOpenMM ky = ry * recipBoxSize[1];
if(ry >= 0) {
for(int n = 0; n < numberOfAtoms; n++) {
for(int n = 0; n < numberOfAtoms; n++)
tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
}
}
else {
for(int n = 0; n < numberOfAtoms; n++) {
for(int n = 0; n < numberOfAtoms; n++)
tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
}
}
for (int rz = lowrz; rz < numRz; rz++) {
if( rz >= 0) {
for( int n = 0; n < numberOfAtoms; n++) {
tab_qxyz[n].real() = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2)).real();
tab_qxyz[n].imag() = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2)).imag();
}
for( int n = 0; n < numberOfAtoms; n++)
tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2));
}
else {
for( int n = 0; n < numberOfAtoms; n++) {
tab_qxyz[n].real() = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2))).real() ;
tab_qxyz[n].imag() = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2))).imag() ;
}
for( int n = 0; n < numberOfAtoms; n++)
tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2)));
}
RealOpenMM cs = 0.0f;
......@@ -400,7 +388,8 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
realSpaceEwaldEnergy = realSpaceEwaldEnergy + atomParameters[atomID1][QIndex]*atomParameters[atomID2][QIndex]*inverseR*erfc(alphaEwald*r);
realSpaceEwaldEnergy =
(RealOpenMM)(realSpaceEwaldEnergy + atomParameters[atomID1][QIndex]*atomParameters[atomID2][QIndex]*inverseR*erfc(alphaEwald*r));
}
}
......@@ -431,9 +420,10 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r;
realSpaceEwaldEnergy = realSpaceEwaldEnergy + atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR);
realSpaceEwaldEnergy =
(RealOpenMM)(realSpaceEwaldEnergy + atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR));
RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
dEdR = dEdR * (erfc(alphaR) + 2.0 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI );
dEdR = (RealOpenMM)(dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
for( int kk = 0; kk < 3; kk++ ){
RealOpenMM force = dEdR*deltaR[0][kk];
......
tests/TestFindExclusions.cpp
View file @ f6352d88
......@@ -109,7 +109,7 @@ int main() {
// This is an exclusion.
ASSERT_EQUAL(0.0, epsilon);
ASSERT(expectedExclusions[particle1].find(particle2) != expectedExclusions[particle2].end());
ASSERT(expectedExclusions[particle1].find(particle2) != expectedExclusions[particle1].end());
}
else {
// This is a 1-4.
......@@ -117,7 +117,7 @@ int main() {
ASSERT_EQUAL_TOL(0.2, chargeProd, 1e-10);
ASSERT_EQUAL_TOL(1.0, sigma, 1e-10);
ASSERT_EQUAL_TOL(0.8, epsilon, 1e-10);
ASSERT(expected14[particle1].find(particle2) != expected14[particle2].end());
ASSERT(expected14[particle1].find(particle2) != expected14[particle1].end());
}
}
}
......
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