Commit f38a47e4 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Cleaned up for Windows compilation

Removed cubic/quartic per-atom parameters from AmoebaHarmonicBondForce
parent 7eaa5d29
...@@ -94,25 +94,23 @@ public: ...@@ -94,25 +94,23 @@ public:
* @param particle1 the index of the first particle connected by the bond * @param particle1 the index of the first particle connected by the bond
* @param particle2 the index of the second particle connected by the bond * @param particle2 the index of the second particle connected by the bond
* @param length the equilibrium length of the bond, measured in nm * @param length the equilibrium length of the bond, measured in nm
* @param cubic k the cubic harmonic force constant for the bond * @param k the quadratic harmonic force constant for the bond
* @param quartic k the quartic harmonic force constant for the bond
* @param quadratic k the quadratic harmonic force constant for the bond
* @return the index of the bond that was added * @return the index of the bond that was added
*/ */
int addBond(int particle1, int particle2, double length, double quadraticK, double cubicK = DEFAULT_GLOBAL_K, double quarticK = DEFAULT_GLOBAL_K );
int addBond(int particle1, int particle2, double length, double quadraticK );
/** /**
* Get the force field parameters for a bond term. * Get the force field parameters for a bond term.
* *
* @param index the index of the bond for which to get parameters * @param index the index of the bond for which to get parameters
* @param particle1 the index of the first particle connected by the bond * @param particle1 the index of the first particle connected by the bond
* @param particle2 the index of the second particle connected by the bond * @param particle2 the index of the second particle connected by the bond
* @param length the equilibrium length of the bond, measured in nm * @param length the equilibrium length of the bond, measured in nm
* @param quadratic k the quadratic harmonic force constant for the bond * @param quadratic k the quadratic harmonic force constant for the bond
* @param cubic k the cubic harmonic force constant for the bond
* @param quartic k the quartic harmonic force constant for the bond
*/ */
void getBondParameters(int index, int& particle1, int& particle2, double& length, double& quadraticK, double& cubicK, double& quarticK ) const;
void getBondParameters(int index, int& particle1, int& particle2, double& length, double& quadraticK ) const;
/** /**
* Set the force field parameters for a bond term. * Set the force field parameters for a bond term.
...@@ -121,18 +119,15 @@ public: ...@@ -121,18 +119,15 @@ public:
* @param particle1 the index of the first particle connected by the bond * @param particle1 the index of the first particle connected by the bond
* @param particle2 the index of the second particle connected by the bond * @param particle2 the index of the second particle connected by the bond
* @param length the equilibrium length of the bond, measured in nm * @param length the equilibrium length of the bond, measured in nm
* @param cubic k the cubic harmonic force constant for the bond * @param k the quadratic harmonic force constant for the bond
* @param quartic k the quartic harmonic force constant for the bond
* @param quadratic k the quadratic harmonic force constant for the bond
*/ */
void setBondParameters(int index, int particle1, int particle2, double length, double quadraticK, double cubicK = DEFAULT_GLOBAL_K, double quarticK = DEFAULT_GLOBAL_K); void setBondParameters(int index, int particle1, int particle2, double length, double quadraticK );
protected: protected:
double _globalQuarticK, _globalCubicK; double _globalQuarticK, _globalCubicK;
ForceImpl* createImpl(); ForceImpl* createImpl();
private: private:
static const double DEFAULT_GLOBAL_K;
class BondInfo; class BondInfo;
// Retarded visual studio compiler complains about being unable to // Retarded visual studio compiler complains about being unable to
...@@ -152,14 +147,13 @@ private: ...@@ -152,14 +147,13 @@ private:
class AmoebaHarmonicBondForce::BondInfo { class AmoebaHarmonicBondForce::BondInfo {
public: public:
int particle1, particle2; int particle1, particle2;
double length, quadraticK, cubicK, quarticK; double length, quadraticK;
BondInfo() { BondInfo() {
particle1 = particle2 = -1; particle1 = particle2 = -1;
length = quadraticK = 0.0; length = quadraticK = 0.0;
cubicK = quarticK = DEFAULT_GLOBAL_K;
} }
BondInfo(int particle1, int particle2, double length, double quadraticK, double cubicK, double quarticK) : BondInfo(int particle1, int particle2, double length, double quadraticK ) :
particle1(particle1), particle2(particle2), length(length), quadraticK(quadraticK), cubicK(cubicK), quarticK(quarticK) { particle1(particle1), particle2(particle2), length(length), quadraticK(quadraticK) {
} }
}; };
......
...@@ -36,33 +36,27 @@ ...@@ -36,33 +36,27 @@
using namespace OpenMM; using namespace OpenMM;
const double AmoebaHarmonicBondForce::DEFAULT_GLOBAL_K = -1234.567891234;
AmoebaHarmonicBondForce::AmoebaHarmonicBondForce() { AmoebaHarmonicBondForce::AmoebaHarmonicBondForce() {
_globalCubicK = _globalQuarticK = 0.0; _globalCubicK = _globalQuarticK = 0.0;
} }
int AmoebaHarmonicBondForce::addBond(int particle1, int particle2, double length, double quadraticK, double cubicK, double quarticK) { int AmoebaHarmonicBondForce::addBond(int particle1, int particle2, double length, double quadraticK) {
bonds.push_back(BondInfo(particle1, particle2, length, quadraticK, cubicK, quarticK)); bonds.push_back(BondInfo(particle1, particle2, length, quadraticK ));
return bonds.size()-1; return bonds.size()-1;
} }
void AmoebaHarmonicBondForce::getBondParameters(int index, int& particle1, int& particle2, double& length, double& quadraticK, double& cubicK, double& quarticK ) const { void AmoebaHarmonicBondForce::getBondParameters(int index, int& particle1, int& particle2, double& length, double& quadraticK ) const {
particle1 = bonds[index].particle1; particle1 = bonds[index].particle1;
particle2 = bonds[index].particle2; particle2 = bonds[index].particle2;
length = bonds[index].length; length = bonds[index].length;
quadraticK = bonds[index].quadraticK; quadraticK = bonds[index].quadraticK;
cubicK = bonds[index].cubicK == DEFAULT_GLOBAL_K ? _globalCubicK : bonds[index].cubicK;
quarticK = bonds[index].quarticK == DEFAULT_GLOBAL_K ? _globalQuarticK : bonds[index].quarticK;
} }
void AmoebaHarmonicBondForce::setBondParameters(int index, int particle1, int particle2, double length, double quadraticK, double cubicK, double quarticK ) { void AmoebaHarmonicBondForce::setBondParameters(int index, int particle1, int particle2, double length, double quadraticK ) {
bonds[index].particle1 = particle1; bonds[index].particle1 = particle1;
bonds[index].particle2 = particle2; bonds[index].particle2 = particle2;
bonds[index].length = length; bonds[index].length = length;
bonds[index].quadraticK = quadraticK; bonds[index].quadraticK = quadraticK;
bonds[index].cubicK = cubicK;
bonds[index].quarticK = quarticK;
} }
void AmoebaHarmonicBondForce::setAmoebaGlobalHarmonicBondCubic(double cubicK ) { void AmoebaHarmonicBondForce::setAmoebaGlobalHarmonicBondCubic(double cubicK ) {
......
...@@ -35,9 +35,6 @@ IF(LOG) ...@@ -35,9 +35,6 @@ IF(LOG)
ENDIF(LOG) ENDIF(LOG)
## ---------------------------------------------------------------------------- ## ----------------------------------------------------------------------------
# message("CUDA_NVCC_FLAGS = ${CUDA_NVCC_FLAGS}")
# INCLUDE(${FINDCUDA_DIR}/FindCuda.cmake)
# message("CUDA_NVCC_FLAGS = ${CUDA_NVCC_FLAGS}")
INCLUDE_DIRECTORIES(${CUDA_INCLUDE}) INCLUDE_DIRECTORIES(${CUDA_INCLUDE})
LINK_DIRECTORIES(${CUDA_TARGET_LINK}) LINK_DIRECTORIES(${CUDA_TARGET_LINK})
FOREACH(subdir ${OPENMM_AMOEBA_SOURCE_SUBDIRS}) FOREACH(subdir ${OPENMM_AMOEBA_SOURCE_SUBDIRS})
...@@ -74,6 +71,6 @@ CUDA_ADD_LIBRARY(${SHARED_TARGET} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES} ...@@ -74,6 +71,6 @@ CUDA_ADD_LIBRARY(${SHARED_TARGET} SHARED ${SOURCE_FILES} ${SOURCE_INCLUDE_FILES}
TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${OPENMM_LIBRARY_NAME}_d optimized ${OPENMM_LIBRARY_NAME} ) TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${OPENMM_LIBRARY_NAME}_d optimized ${OPENMM_LIBRARY_NAME} )
TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${OPENMM_LIBRARY_NAME}Cuda_d optimized ${OPENMM_LIBRARY_NAME}Cuda ) TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${OPENMM_LIBRARY_NAME}Cuda_d optimized ${OPENMM_LIBRARY_NAME}Cuda )
TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${OPENMM_AMOEBA_LIBRARY_NAME}_d optimized ${OPENMM_AMOEBA_LIBRARY_NAME} ) TARGET_LINK_LIBRARIES(${SHARED_TARGET} debug ${OPENMM_AMOEBA_LIBRARY_NAME}_d optimized ${OPENMM_AMOEBA_LIBRARY_NAME} )
SET_TARGET_PROPERTIES(${SHARED_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMMCUDAAMOEBA_BUILDING_SHARED_LIBRARY") SET_TARGET_PROPERTIES(${SHARED_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMMCUDA_BUILDING_SHARED_LIBRARY -DOPENMMCUDAAMOEBA_BUILDING_SHARED_LIBRARY")
INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${SHARED_TARGET}) INSTALL_TARGETS(/lib/plugins RUNTIME_DIRECTORY /lib/plugins ${SHARED_TARGET})
...@@ -96,7 +96,7 @@ public: ...@@ -96,7 +96,7 @@ public:
* *
* @return amoebaGpuContext * @return amoebaGpuContext
*/ */
amoebaGpuContext getAmoebaGpu( void ) const; amoebaGpuContext OPENMMCUDA_EXPORT getAmoebaGpu( void ) const;
/** /**
* Set accessor for LocalForcesKernel * Set accessor for LocalForcesKernel
......
...@@ -73,23 +73,21 @@ void CudaCalcAmoebaHarmonicBondForceKernel::initialize(const System& system, con ...@@ -73,23 +73,21 @@ void CudaCalcAmoebaHarmonicBondForceKernel::initialize(const System& system, con
std::vector<int> particle2(numBonds); std::vector<int> particle2(numBonds);
std::vector<float> length(numBonds); std::vector<float> length(numBonds);
std::vector<float> quadratic(numBonds); std::vector<float> quadratic(numBonds);
std::vector<float> cubic(numBonds);
std::vector<float> quartic(numBonds);
for (int i = 0; i < numBonds; i++) { for (int i = 0; i < numBonds; i++) {
int particle1Index, particle2Index; int particle1Index, particle2Index;
double lengthValue, kValue, cubicValue, quarticValue; double lengthValue, kValue;
force.getBondParameters(i, particle1Index, particle2Index, lengthValue, kValue, cubicValue, quarticValue); force.getBondParameters(i, particle1Index, particle2Index, lengthValue, kValue );
particle1[i] = particle1Index; particle1[i] = particle1Index;
particle2[i] = particle2Index; particle2[i] = particle2Index;
length[i] = static_cast<float>( lengthValue ); length[i] = static_cast<float>( lengthValue );
quadratic[i] = static_cast<float>( kValue ); quadratic[i] = static_cast<float>( kValue );
cubic[i] = static_cast<float>( cubicValue );
quartic[i] = static_cast<float>( quarticValue );
} }
gpuSetAmoebaBondParameters( data.getAmoebaGpu(), particle1, particle2, length, quadratic, cubic, quartic); gpuSetAmoebaBondParameters( data.getAmoebaGpu(), particle1, particle2, length, quadratic,
static_cast<float>(force.getAmoebaGlobalHarmonicBondCubic()),
static_cast<float>(force.getAmoebaGlobalHarmonicBondQuartic()) );
} }
void CudaCalcAmoebaHarmonicBondForceKernel::executeForces(ContextImpl& context) { void CudaCalcAmoebaHarmonicBondForceKernel::executeForces(ContextImpl& context) {
...@@ -132,10 +130,10 @@ void CudaCalcAmoebaHarmonicAngleForceKernel::initialize(const System& system, co ...@@ -132,10 +130,10 @@ void CudaCalcAmoebaHarmonicAngleForceKernel::initialize(const System& system, co
k[i] = static_cast<float>( kQuadratic ); k[i] = static_cast<float>( kQuadratic );
} }
gpuSetAmoebaAngleParameters(data.getAmoebaGpu(), particle1, particle2, particle3, angle, k, gpuSetAmoebaAngleParameters(data.getAmoebaGpu(), particle1, particle2, particle3, angle, k,
force.getAmoebaGlobalHarmonicAngleCubic(), static_cast<float>(force.getAmoebaGlobalHarmonicAngleCubic()),
force.getAmoebaGlobalHarmonicAngleQuartic(), static_cast<float>(force.getAmoebaGlobalHarmonicAngleQuartic()),
force.getAmoebaGlobalHarmonicAnglePentic(), static_cast<float>(force.getAmoebaGlobalHarmonicAnglePentic()),
force.getAmoebaGlobalHarmonicAngleSextic() ); static_cast<float>(force.getAmoebaGlobalHarmonicAngleSextic()) );
} }
...@@ -182,10 +180,10 @@ void CudaCalcAmoebaHarmonicInPlaneAngleForceKernel::initialize(const System& sys ...@@ -182,10 +180,10 @@ void CudaCalcAmoebaHarmonicInPlaneAngleForceKernel::initialize(const System& sys
k[i] = static_cast<float>( kQuadratic ); k[i] = static_cast<float>( kQuadratic );
} }
gpuSetAmoebaInPlaneAngleParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, angle, k, gpuSetAmoebaInPlaneAngleParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, angle, k,
force.getAmoebaGlobalHarmonicInPlaneAngleCubic(), static_cast<float>( force.getAmoebaGlobalHarmonicInPlaneAngleCubic()),
force.getAmoebaGlobalHarmonicInPlaneAngleQuartic(), static_cast<float>( force.getAmoebaGlobalHarmonicInPlaneAngleQuartic()),
force.getAmoebaGlobalHarmonicInPlaneAnglePentic(), static_cast<float>( force.getAmoebaGlobalHarmonicInPlaneAnglePentic()),
force.getAmoebaGlobalHarmonicInPlaneAngleSextic() ); static_cast<float>( force.getAmoebaGlobalHarmonicInPlaneAngleSextic() ) );
} }
...@@ -237,9 +235,9 @@ void CudaCalcAmoebaTorsionForceKernel::initialize(const System& system, const Am ...@@ -237,9 +235,9 @@ void CudaCalcAmoebaTorsionForceKernel::initialize(const System& system, const Am
force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], torsionParameter1, torsionParameter2, torsionParameter3 ); force.getTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], torsionParameter1, torsionParameter2, torsionParameter3 );
for ( unsigned int jj = 0; jj < 3; jj++) { for ( unsigned int jj = 0; jj < 3; jj++) {
torsionParameters1F[jj] = torsionParameter1[jj]; torsionParameters1F[jj] = static_cast<float>(torsionParameter1[jj]);
torsionParameters2F[jj] = torsionParameter2[jj]; torsionParameters2F[jj] = static_cast<float>(torsionParameter2[jj]);
torsionParameters3F[jj] = torsionParameter3[jj]; torsionParameters3F[jj] = static_cast<float>(torsionParameter3[jj]);
} }
torsionParameters1[i] = torsionParameters1F; torsionParameters1[i] = torsionParameters1F;
torsionParameters2[i] = torsionParameters2F; torsionParameters2[i] = torsionParameters2F;
...@@ -290,7 +288,7 @@ void CudaCalcAmoebaPiTorsionForceKernel::initialize(const System& system, const ...@@ -290,7 +288,7 @@ void CudaCalcAmoebaPiTorsionForceKernel::initialize(const System& system, const
double torsionKParameter; double torsionKParameter;
force.getPiTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], particle5[i], particle6[i], torsionKParameter); force.getPiTorsionParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], particle5[i], particle6[i], torsionKParameter);
torsionKParameters[i] = torsionKParameter; torsionKParameters[i] = static_cast<float>(torsionKParameter);
} }
gpuSetAmoebaPiTorsionParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, particle5, particle6, torsionKParameters); gpuSetAmoebaPiTorsionParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, particle5, particle6, torsionKParameters);
} }
...@@ -335,10 +333,10 @@ void CudaCalcAmoebaStretchBendForceKernel::initialize(const System& system, cons ...@@ -335,10 +333,10 @@ void CudaCalcAmoebaStretchBendForceKernel::initialize(const System& system, cons
double lengthAB, lengthCB, angle, k; double lengthAB, lengthCB, angle, k;
force.getStretchBendParameters(i, particle1[i], particle2[i], particle3[i], lengthAB, lengthCB, angle, k); force.getStretchBendParameters(i, particle1[i], particle2[i], particle3[i], lengthAB, lengthCB, angle, k);
lengthABParameters[i] = lengthAB; lengthABParameters[i] = static_cast<float>(lengthAB);
lengthCBParameters[i] = lengthCB; lengthCBParameters[i] = static_cast<float>(lengthCB);
angleParameters[i] = angle; angleParameters[i] = static_cast<float>(angle);
kParameters[i] = k; kParameters[i] = static_cast<float>(k);
} }
gpuSetAmoebaStretchBendParameters(data.getAmoebaGpu(), particle1, particle2, particle3, lengthABParameters, lengthCBParameters, angleParameters, kParameters); gpuSetAmoebaStretchBendParameters(data.getAmoebaGpu(), particle1, particle2, particle3, lengthABParameters, lengthCBParameters, angleParameters, kParameters);
...@@ -381,13 +379,13 @@ void CudaCalcAmoebaOutOfPlaneBendForceKernel::initialize(const System& system, c ...@@ -381,13 +379,13 @@ void CudaCalcAmoebaOutOfPlaneBendForceKernel::initialize(const System& system, c
double k; double k;
force.getOutOfPlaneBendParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], k); force.getOutOfPlaneBendParameters(i, particle1[i], particle2[i], particle3[i], particle4[i], k);
kParameters[i] = k; kParameters[i] = static_cast<float>(k);
} }
gpuSetAmoebaOutOfPlaneBendParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, kParameters, gpuSetAmoebaOutOfPlaneBendParameters(data.getAmoebaGpu(), particle1, particle2, particle3, particle4, kParameters,
force.getAmoebaGlobalOutOfPlaneBendCubic(), static_cast<float>( force.getAmoebaGlobalOutOfPlaneBendCubic()),
force.getAmoebaGlobalOutOfPlaneBendQuartic(), static_cast<float>( force.getAmoebaGlobalOutOfPlaneBendQuartic()),
force.getAmoebaGlobalOutOfPlaneBendPentic(), static_cast<float>( force.getAmoebaGlobalOutOfPlaneBendPentic()),
force.getAmoebaGlobalOutOfPlaneBendSextic() ); static_cast<float>( force.getAmoebaGlobalOutOfPlaneBendSextic() ) );
} }
......
...@@ -29,16 +29,18 @@ ...@@ -29,16 +29,18 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. * * USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
//using namespace std;
#include "cudaKernels.h" #include "cudaKernels.h"
#include "amoebaCudaKernels.h" #include "amoebaCudaKernels.h"
// for some reason, these are not being included w/ cudaKernels.h on Windows
extern void OPENMMCUDA_EXPORT SetCalculateObcGbsaForces2Sim(gpuContext gpu);
extern void OPENMMCUDA_EXPORT SetForcesSim(gpuContext gpu);
#include <sstream> #include <sstream>
#include <limits> #include <limits>
#ifdef WIN32 #ifdef WIN32
#include <windows.h> // #include <windows.h>
#endif #endif
#define DUMP_PARAMETERS 0 #define DUMP_PARAMETERS 0
...@@ -207,9 +209,11 @@ void gpuPrintCudaAmoebaGmxSimulation(amoebaGpuContext amoebaGpu, FILE* log ) ...@@ -207,9 +209,11 @@ void gpuPrintCudaAmoebaGmxSimulation(amoebaGpuContext amoebaGpu, FILE* log )
if( amoebaGpu->psAmoebaBondParameter)(void) fprintf( log, "\n" ); if( amoebaGpu->psAmoebaBondParameter)(void) fprintf( log, "\n" );
gpuPrintCudaStreamInt4( amoebaGpu->psAmoebaBondID, log ); gpuPrintCudaStreamInt4( amoebaGpu->psAmoebaBondID, log );
gpuPrintCudaStreamFloat4( amoebaGpu->psAmoebaBondParameter, log ); gpuPrintCudaStreamFloat2( amoebaGpu->psAmoebaBondParameter, log );
(void) fprintf( log, " amoebaBonds %u\n", amoebaGpu->amoebaSim.amoebaBonds ); (void) fprintf( log, " amoebaBonds %u\n", amoebaGpu->amoebaSim.amoebaBonds );
(void) fprintf( log, " amoebaBond_offset %u\n", amoebaGpu->amoebaSim.amoebaBond_offset ); (void) fprintf( log, " amoebaBond_offset %u\n", amoebaGpu->amoebaSim.amoebaBond_offset );
(void) fprintf( log, " cubic %14.7e\n", amoebaGpu->amoebaSim.amoebaBondCubicParameter);
(void) fprintf( log, " quartic %14.7e\n", amoebaGpu->amoebaSim.amoebaBondQuarticicParameter);
(void) fprintf( log, " pAmoebaBondID %p\n", amoebaGpu->amoebaSim.pAmoebaBondID ); (void) fprintf( log, " pAmoebaBondID %p\n", amoebaGpu->amoebaSim.pAmoebaBondID );
(void) fprintf( log, " pAmoebaBondParameter %p\n", amoebaGpu->amoebaSim.pAmoebaBondParameter ); (void) fprintf( log, " pAmoebaBondParameter %p\n", amoebaGpu->amoebaSim.pAmoebaBondParameter );
...@@ -415,18 +419,22 @@ void gpuPrintCudaAmoebaGmxSimulation(amoebaGpuContext amoebaGpu, FILE* log ) ...@@ -415,18 +419,22 @@ void gpuPrintCudaAmoebaGmxSimulation(amoebaGpuContext amoebaGpu, FILE* log )
extern "C" extern "C"
void gpuSetAmoebaBondParameters(amoebaGpuContext amoebaGpu, const std::vector<int>& particles1, const std::vector<int>& particles2, void gpuSetAmoebaBondParameters(amoebaGpuContext amoebaGpu, const std::vector<int>& particles1, const std::vector<int>& particles2,
const std::vector<float>& length, const std::vector<float>& k, const std::vector<float>& cubic, const std::vector<float>& length, const std::vector<float>& k, float cubic, float quartic)
const std::vector<float>& quartic)
{ {
_gpuContext* gpu = amoebaGpu->gpuContext; _gpuContext* gpu = amoebaGpu->gpuContext;
int bonds = particles1.size(); int bonds = particles1.size();
amoebaGpu->amoebaSim.amoebaBonds = bonds; amoebaGpu->amoebaSim.amoebaBonds = bonds;
CUDAStream<int4>* psBondID = new CUDAStream<int4>(bonds, 1, "AmoebaBondID");
amoebaGpu->psAmoebaBondID = psBondID; CUDAStream<int4>* psBondID = new CUDAStream<int4>(bonds, 1, "AmoebaBondID");
amoebaGpu->amoebaSim.pAmoebaBondID = psBondID->_pDevStream[0]; amoebaGpu->psAmoebaBondID = psBondID;
CUDAStream<float4>* psBondParameter = new CUDAStream<float4>(bonds, 1, "AmoebaBondParameter"); amoebaGpu->amoebaSim.pAmoebaBondID = psBondID->_pDevStream[0];
amoebaGpu->psAmoebaBondParameter = psBondParameter;
amoebaGpu->amoebaSim.pAmoebaBondParameter = psBondParameter->_pDevStream[0]; CUDAStream<float2>* psBondParameter = new CUDAStream<float2>(bonds, 1, "AmoebaBondParameter");
amoebaGpu->psAmoebaBondParameter = psBondParameter;
amoebaGpu->amoebaSim.pAmoebaBondParameter = psBondParameter->_pDevStream[0];
amoebaGpu->amoebaSim.amoebaBondCubicParameter = cubic;
amoebaGpu->amoebaSim.amoebaBondQuarticicParameter = quartic;
for (int i = 0; i < bonds; i++) for (int i = 0; i < bonds; i++)
{ {
(*psBondID)[i].x = particles1[i]; (*psBondID)[i].x = particles1[i];
...@@ -435,8 +443,6 @@ void gpuSetAmoebaBondParameters(amoebaGpuContext amoebaGpu, const std::vector<in ...@@ -435,8 +443,6 @@ void gpuSetAmoebaBondParameters(amoebaGpuContext amoebaGpu, const std::vector<in
(*psBondID)[i].w = gpu->pOutputBufferCounter[(*psBondID)[i].y]++; (*psBondID)[i].w = gpu->pOutputBufferCounter[(*psBondID)[i].y]++;
(*psBondParameter)[i].x = length[i]; (*psBondParameter)[i].x = length[i];
(*psBondParameter)[i].y = k[i]; (*psBondParameter)[i].y = k[i];
(*psBondParameter)[i].z = cubic[i];
(*psBondParameter)[i].w = quartic[i];
#undef DUMP_PARAMETERS #undef DUMP_PARAMETERS
#define DUMP_PARAMETERS 5 #define DUMP_PARAMETERS 5
...@@ -444,7 +450,7 @@ void gpuSetAmoebaBondParameters(amoebaGpuContext amoebaGpu, const std::vector<in ...@@ -444,7 +450,7 @@ void gpuSetAmoebaBondParameters(amoebaGpuContext amoebaGpu, const std::vector<in
if( amoebaGpu->log && (i < DUMP_PARAMETERS || i > bonds - (DUMP_PARAMETERS + 1) ) ) if( amoebaGpu->log && (i < DUMP_PARAMETERS || i > bonds - (DUMP_PARAMETERS + 1) ) )
fprintf( amoebaGpu->log, "Bonds: %5d [%5d %5d %5d %5d] L=%14.7e k[%14.7e %14.7e %14.7e] [%5d %5d]\n", fprintf( amoebaGpu->log, "Bonds: %5d [%5d %5d %5d %5d] L=%14.7e k[%14.7e %14.7e %14.7e] [%5d %5d]\n",
i, (*psBondID)[i].x, (*psBondID)[i].y, (*psBondID)[i].z, (*psBondID)[i].w, i, (*psBondID)[i].x, (*psBondID)[i].y, (*psBondID)[i].z, (*psBondID)[i].w,
(*psBondParameter)[i].x, (*psBondParameter)[i].y, (*psBondParameter)[i].z, (*psBondParameter)[i].w, (*psBondParameter)[i].x, (*psBondParameter)[i].y, cubic, quartic,
gpu->pOutputBufferCounter[(*psBondID)[i].x], gpu->pOutputBufferCounter[(*psBondID)[i].x],
gpu->pOutputBufferCounter[(*psBondID)[i].y] ); gpu->pOutputBufferCounter[(*psBondID)[i].y] );
#endif #endif
...@@ -938,7 +944,7 @@ void gpuSetAmoebaTorsionTorsionGrids(amoebaGpuContext amoebaGpu, const std::vect ...@@ -938,7 +944,7 @@ void gpuSetAmoebaTorsionTorsionGrids(amoebaGpuContext amoebaGpu, const std::vect
(void) fprintf( amoebaGpu->log, "totalGridEntries=%u totalFloat4 entries=%u\n", totalGridEntries, totalEntries ); (void) fprintf( amoebaGpu->log, "totalGridEntries=%u totalFloat4 entries=%u\n", totalGridEntries, totalEntries );
} }
int index = 0; unsigned int index = 0;
for (unsigned int ii = 0; ii < floatGrids.size(); ii++) { for (unsigned int ii = 0; ii < floatGrids.size(); ii++) {
for (unsigned int jj = 0; jj < floatGrids[ii].size(); jj++) { for (unsigned int jj = 0; jj < floatGrids[ii].size(); jj++) {
for (unsigned int kk = 0; kk < floatGrids[ii][jj].size(); kk++) { for (unsigned int kk = 0; kk < floatGrids[ii][jj].size(); kk++) {
...@@ -1053,7 +1059,9 @@ void gpuSetAmoebaBondOffsets(amoebaGpuContext amoebaGpu ) ...@@ -1053,7 +1059,9 @@ void gpuSetAmoebaBondOffsets(amoebaGpuContext amoebaGpu )
amoebaGpu->amoebaSim.amoebaTorsionTorsion_offset = amoebaGpu->amoebaSim.amoebaOutOfPlaneBend_offset + amoebaGpu->amoebaSim.amoebaTorsionTorsion_offset = amoebaGpu->amoebaSim.amoebaOutOfPlaneBend_offset +
(amoebaGpu->psAmoebaTorsionTorsionID1 ? amoebaGpu->psAmoebaTorsionTorsionID1->_stride : 0); (amoebaGpu->psAmoebaTorsionTorsionID1 ? amoebaGpu->psAmoebaTorsionTorsionID1->_stride : 0);
gpu->sim.localForces_threads_per_block = (std::max(amoebaGpu->amoebaSim.amoebaTorsionTorsion_offset, gpu->sim.customBonds) / gpu->sim.blocks + 15) & 0xfffffff0; //gpu->sim.localForces_threads_per_block = (std::max(amoebaGpu->amoebaSim.amoebaTorsionTorsion_offset, gpu->sim.customBonds) / gpu->sim.blocks + 15) & 0xfffffff0;
unsigned int maxI = (amoebaGpu->amoebaSim.amoebaTorsionTorsion_offset > gpu->sim.customBonds) ? amoebaGpu->amoebaSim.amoebaTorsionTorsion_offset : gpu->sim.customBonds;
gpu->sim.localForces_threads_per_block = (maxI/gpu->sim.blocks + 15) & 0xfffffff0;
if (gpu->sim.localForces_threads_per_block > gpu->sim.max_localForces_threads_per_block) if (gpu->sim.localForces_threads_per_block > gpu->sim.max_localForces_threads_per_block)
gpu->sim.localForces_threads_per_block = gpu->sim.max_localForces_threads_per_block; gpu->sim.localForces_threads_per_block = gpu->sim.max_localForces_threads_per_block;
if (gpu->sim.localForces_threads_per_block < 1) if (gpu->sim.localForces_threads_per_block < 1)
...@@ -1495,7 +1503,7 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect ...@@ -1495,7 +1503,7 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
amoebaGpu->amoebaSim.dielec = 1.0f; amoebaGpu->amoebaSim.dielec = 1.0f;
} }
for( unsigned int ii = 0; ii < charges.size(); ii++ ){ for( int ii = 0; ii < static_cast<int>(charges.size()); ii++ ){
// axis type & multipole particles ids // axis type & multipole particles ids
...@@ -2140,16 +2148,16 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec ...@@ -2140,16 +2148,16 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec
std::vector<int> minCellIndex( dim + 1 ); std::vector<int> minCellIndex( dim + 1 );
std::vector<int> maxCellIndex( dim + 1 ); std::vector<int> maxCellIndex( dim + 1 );
for (int ii = 0; ii <= dim; ii++) for(unsigned int ii = 0; ii <= dim; ii++)
{ {
minCellIndex[ii] = paddedAtoms + 1; minCellIndex[ii] = paddedAtoms + 1;
maxCellIndex[ii] = 0; maxCellIndex[ii] = 0;
} }
for (unsigned int atom1 = 0; atom1 < actualAtoms; atom1++) for(unsigned int atom1 = 0; atom1 < actualAtoms; atom1++)
{ {
int x = atom1/grid; int x = atom1/grid;
for (int jj = 0; jj < exclusions[atom1].size(); jj++ ) for ( unsigned int jj = 0; jj < exclusions[atom1].size(); jj++ )
{ {
if( exclusions[atom1][jj] > maxCellIndex[x] ) if( exclusions[atom1][jj] > maxCellIndex[x] )
{ {
...@@ -2166,7 +2174,7 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec ...@@ -2166,7 +2174,7 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec
if( debugOn && amoebaGpu->log ){ if( debugOn && amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s min/max cell indices:\n", methodName.c_str() ); (void) fprintf( amoebaGpu->log, "%s min/max cell indices:\n", methodName.c_str() );
for (int ii = 0; ii < dim; ii++) for ( unsigned int ii = 0; ii < dim; ii++)
{ {
(void) fprintf( amoebaGpu->log, "%6d [%6d %6d]\n", ii, minCellIndex[ii], maxCellIndex[ii] ); (void) fprintf( amoebaGpu->log, "%6d [%6d %6d]\n", ii, minCellIndex[ii], maxCellIndex[ii] );
} }
...@@ -2180,14 +2188,14 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec ...@@ -2180,14 +2188,14 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec
amoebaGpu->amoebaSim.pVdwExclusionIndicesIndex = psVdwExclusionIndicesIndex->_pDevStream[0]; amoebaGpu->amoebaSim.pVdwExclusionIndicesIndex = psVdwExclusionIndicesIndex->_pDevStream[0];
//memset( amoebaGpu->psVdwExclusionIndicesIndex->_pSysStream[0], 0, sizeof(cells)*sizeof(int) ); //memset( amoebaGpu->psVdwExclusionIndicesIndex->_pSysStream[0], 0, sizeof(cells)*sizeof(int) );
for (int ii = 0; ii < cells; ii++) for (unsigned int ii = 0; ii < cells; ii++)
{ {
amoebaGpu->psVdwExclusionIndicesIndex->_pSysStream[0][ii] = -1; amoebaGpu->psVdwExclusionIndicesIndex->_pSysStream[0][ii] = -1;
} }
int numWithExclusionIndices = 0; int numWithExclusionIndices = 0;
int gridOffset = grid - 1; int gridOffset = grid - 1;
int lastBlock = (paddedAtoms > amoebaGpu->gpuContext->natoms) ? (amoebaGpu->gpuContext->natoms)/grid : -1; int lastBlock = (static_cast<int>(paddedAtoms) > amoebaGpu->gpuContext->natoms) ? (amoebaGpu->gpuContext->natoms)/grid : -1;
for (int ii = 0; ii < cells; ii++) for (unsigned int ii = 0; ii < cells; ii++)
{ {
unsigned int x, y, exclusion; unsigned int x, y, exclusion;
decodeCell( pWorkList[ii], &x, &y, &exclusion ); decodeCell( pWorkList[ii], &x, &y, &exclusion );
...@@ -2201,7 +2209,7 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec ...@@ -2201,7 +2209,7 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec
// ii, x, y, xAtomMin, xAtomMax, minCellIndex[y], maxCellIndex[y], numWithExclusionIndices, lastBlock ); // ii, x, y, xAtomMin, xAtomMax, minCellIndex[y], maxCellIndex[y], numWithExclusionIndices, lastBlock );
} }
} }
for (int ii = 0; ii < cells; ii++) for ( unsigned int ii = 0; ii < cells; ii++)
{ {
if( amoebaGpu->psVdwExclusionIndicesIndex->_pSysStream[0][ii] == -1 ) if( amoebaGpu->psVdwExclusionIndicesIndex->_pSysStream[0][ii] == -1 )
{ {
...@@ -2213,7 +2221,7 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec ...@@ -2213,7 +2221,7 @@ void amoebaGpuBuildVdwExclusionList( amoebaGpuContext amoebaGpu, const std::vec
if( debugOn && amoebaGpu->log ){ if( debugOn && amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s %d cells w/ exclusions\n", methodName.c_str(), numWithExclusionIndices ); (void) fprintf( amoebaGpu->log, "%s %d cells w/ exclusions\n", methodName.c_str(), numWithExclusionIndices );
for (int ii = 0; ii < cells; ii++) for (unsigned int ii = 0; ii < cells; ii++)
{ {
unsigned int x, y, exclusion; unsigned int x, y, exclusion;
decodeCell( pWorkList[ii], &x, &y, &exclusion ); decodeCell( pWorkList[ii], &x, &y, &exclusion );
...@@ -2774,10 +2782,10 @@ static int matchMaps( std::string idString, MapIntFloat* map1, MapIntFloat* map ...@@ -2774,10 +2782,10 @@ static int matchMaps( std::string idString, MapIntFloat* map1, MapIntFloat* map
static void getScalingDegrees( amoebaGpuContext amoebaGpu, unsigned int particleI, unsigned int particleJ, int* covalentDegree, int* polarizationDegree ) static void getScalingDegrees( amoebaGpuContext amoebaGpu, unsigned int particleI, unsigned int particleJ, int* covalentDegree, int* polarizationDegree )
{ {
int particlesOffset = particleI*amoebaGpu->maxCovalentDegreeSz; int particlesOffset = particleI*amoebaGpu->maxCovalentDegreeSz;
int minCovalentIndex = amoebaGpu->psCovalentDegree->_pSysStream[0][particlesOffset]; unsigned int minCovalentIndex = static_cast<unsigned int>(amoebaGpu->psCovalentDegree->_pSysStream[0][particlesOffset]);
int minCovalentPolarizationIndex = amoebaGpu->psPolarizationDegree->_pSysStream[0][particlesOffset]; unsigned int minCovalentPolarizationIndex = static_cast<unsigned int>( amoebaGpu->psPolarizationDegree->_pSysStream[0][particlesOffset]);
if( particleJ < minCovalentIndex || particleJ > (minCovalentIndex + amoebaGpu->maxCovalentDegreeSz) ){ if( particleJ < minCovalentIndex || particleJ > (minCovalentIndex + amoebaGpu->maxCovalentDegreeSz) ){
*covalentDegree = 0; *covalentDegree = 0;
...@@ -2923,13 +2931,13 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu ) ...@@ -2923,13 +2931,13 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu )
minCellIndex.resize( dim + 1 ); minCellIndex.resize( dim + 1 );
maxCellIndex.resize( dim + 1 ); maxCellIndex.resize( dim + 1 );
for (int ii = 0; ii <= dim; ii++) for ( unsigned int ii = 0; ii <= dim; ii++)
{ {
minCellIndex[ii] = paddedAtoms + 1; minCellIndex[ii] = paddedAtoms + 1;
maxCellIndex[ii] = 0; maxCellIndex[ii] = 0;
} }
for (int atom1 = 0; atom1 < paddedAtoms; atom1++) for (unsigned int atom1 = 0; atom1 < paddedAtoms; atom1++)
{ {
int x = atom1/grid; int x = atom1/grid;
int particlesOffset = atom1*amoebaGpu->maxCovalentDegreeSz; int particlesOffset = atom1*amoebaGpu->maxCovalentDegreeSz;
...@@ -2962,7 +2970,7 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu ) ...@@ -2962,7 +2970,7 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu )
if( debugOn && amoebaGpu->log ){ if( debugOn && amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s min/max cell indices:\n", methodName.c_str() ); (void) fprintf( amoebaGpu->log, "%s min/max cell indices:\n", methodName.c_str() );
for (int ii = 0; ii < dim; ii++) for (unsigned int ii = 0; ii < dim; ii++)
{ {
(void) fprintf( amoebaGpu->log, "%6d [%6d %6d]\n", ii, minCellIndex[ii], maxCellIndex[ii] ); (void) fprintf( amoebaGpu->log, "%6d [%6d %6d]\n", ii, minCellIndex[ii], maxCellIndex[ii] );
} }
...@@ -2978,8 +2986,8 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu ) ...@@ -2978,8 +2986,8 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu )
memset( amoebaGpu->psScalingIndicesIndex->_pSysStream[0], 0, sizeof( cells)*sizeof( unsigned int) ); memset( amoebaGpu->psScalingIndicesIndex->_pSysStream[0], 0, sizeof( cells)*sizeof( unsigned int) );
int numWithScalingIndices = 0; int numWithScalingIndices = 0;
int gridOffset = grid - 1; int gridOffset = grid - 1;
int lastBlock = (paddedAtoms > amoebaGpu->gpuContext->natoms) ? (amoebaGpu->gpuContext->natoms)/grid : -1; int lastBlock = (static_cast<int>(paddedAtoms) > amoebaGpu->gpuContext->natoms) ? (amoebaGpu->gpuContext->natoms)/grid : -1;
for (int ii = 0; ii < cells; ii++) for (unsigned int ii = 0; ii < cells; ii++)
{ {
unsigned int x, y, exclusion; unsigned int x, y, exclusion;
decodeCell( pWorkList[ii], &x, &y, &exclusion ); decodeCell( pWorkList[ii], &x, &y, &exclusion );
...@@ -3019,7 +3027,7 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu ) ...@@ -3019,7 +3027,7 @@ void amoebaGpuBuildScalingList( amoebaGpuContext amoebaGpu )
if( debugOn && amoebaGpu->log ){ if( debugOn && amoebaGpu->log ){
(void) fprintf( amoebaGpu->log, "%s %d cells w/ exclusions\n", (void) fprintf( amoebaGpu->log, "%s %d cells w/ exclusions\n",
methodName.c_str(), numWithScalingIndices ); methodName.c_str(), numWithScalingIndices );
for (int ii = 0; ii < cells; ii++) for (unsigned int ii = 0; ii < cells; ii++)
{ {
unsigned int x, y, exclusion; unsigned int x, y, exclusion;
decodeCell( pWorkList[ii], &x, &y, &exclusion ); decodeCell( pWorkList[ii], &x, &y, &exclusion );
...@@ -3703,7 +3711,7 @@ void cudaWriteVectorOfDoubleVectorsToFile( char* fname, std::vector<int>& fileId ...@@ -3703,7 +3711,7 @@ void cudaWriteVectorOfDoubleVectorsToFile( char* fname, std::vector<int>& fileId
for ( unsigned int ii = 0; ii < outputVector.size(); ii++ ){ for ( unsigned int ii = 0; ii < outputVector.size(); ii++ ){
int index = 0; int index = 0;
for ( unsigned int jj = 0; jj < outputVector[ii].size(); jj++ ){ for ( unsigned int jj = 0; jj < outputVector[ii].size(); jj++ ){
values[index++] = outputVector[ii][jj]; values[index++] = static_cast<float>(outputVector[ii][jj]);
} }
printValues( filePtr, static_cast<int>(ii), index, values ); printValues( filePtr, static_cast<int>(ii), index, values );
} }
......
...@@ -46,7 +46,9 @@ struct cudaAmoebaGmxSimulation { ...@@ -46,7 +46,9 @@ struct cudaAmoebaGmxSimulation {
unsigned int amoebaBonds; // Number of bonds unsigned int amoebaBonds; // Number of bonds
int4* pAmoebaBondID; // Bond atom and output buffer IDs int4* pAmoebaBondID; // Bond atom and output buffer IDs
float4* pAmoebaBondParameter; // Bond parameters float2* pAmoebaBondParameter; // Bond parameters
float amoebaBondCubicParameter; // cubic bond parameters
float amoebaBondQuarticicParameter; // quartic bond parameters
unsigned int amoebaBond_offset; // Offset to end of bonds unsigned int amoebaBond_offset; // Offset to end of bonds
unsigned int amoebaAngles; // Number of bond angles unsigned int amoebaAngles; // Number of bond angles
......
...@@ -88,7 +88,7 @@ struct _amoebaGpuContext { ...@@ -88,7 +88,7 @@ struct _amoebaGpuContext {
int maxCovalentDegreeSz; int maxCovalentDegreeSz;
CUDAStream<int4>* psAmoebaBondID; CUDAStream<int4>* psAmoebaBondID;
CUDAStream<float4>* psAmoebaBondParameter; CUDAStream<float2>* psAmoebaBondParameter;
CUDAStream<int4>* psAmoebaAngleID1; CUDAStream<int4>* psAmoebaAngleID1;
CUDAStream<int2>* psAmoebaAngleID2; CUDAStream<int2>* psAmoebaAngleID2;
...@@ -254,7 +254,7 @@ void amoebaGpuBuildScalingList( amoebaGpuContext gpu ); ...@@ -254,7 +254,7 @@ void amoebaGpuBuildScalingList( amoebaGpuContext gpu );
extern "C" extern "C"
void gpuSetAmoebaBondParameters(amoebaGpuContext gpu, const std::vector<int>& atom1, const std::vector<int>& atom2, void gpuSetAmoebaBondParameters(amoebaGpuContext gpu, const std::vector<int>& atom1, const std::vector<int>& atom2,
const std::vector<float>& length, const std::vector<float>& k, const std::vector<float>& cubic, const std::vector<float>& quartic); const std::vector<float>& length, const std::vector<float>& k, float cubic, float quartic);
extern "C" extern "C"
void gpuSetAmoebaAngleParameters(amoebaGpuContext gpu, const std::vector<int>& atom1, const std::vector<int>& atom2, const std::vector<int>& atom3, void gpuSetAmoebaAngleParameters(amoebaGpuContext gpu, const std::vector<int>& atom1, const std::vector<int>& atom2, const std::vector<int>& atom3,
......
...@@ -317,7 +317,7 @@ void kCalculateAmoebaLocalForces_kernel() ...@@ -317,7 +317,7 @@ void kCalculateAmoebaLocalForces_kernel()
int4 atom = cAmoebaSim.pAmoebaBondID[pos]; int4 atom = cAmoebaSim.pAmoebaBondID[pos];
float4 atomA = cSim.pPosq[atom.x]; float4 atomA = cSim.pPosq[atom.x];
float4 atomB = cSim.pPosq[atom.y]; float4 atomB = cSim.pPosq[atom.y];
float4 bond = cAmoebaSim.pAmoebaBondParameter[pos]; float2 bond = cAmoebaSim.pAmoebaBondParameter[pos];
float dx = atomB.x - atomA.x; float dx = atomB.x - atomA.x;
float dy = atomB.y - atomA.y; float dy = atomB.y - atomA.y;
float dz = atomB.z - atomA.z; float dz = atomB.z - atomA.z;
...@@ -330,8 +330,12 @@ void kCalculateAmoebaLocalForces_kernel() ...@@ -330,8 +330,12 @@ void kCalculateAmoebaLocalForces_kernel()
dEdR = (r > 0.0f) ? (dEdR / r) : 0.0f; dEdR = (r > 0.0f) ? (dEdR / r) : 0.0f;
#else #else
float deltaIdeal2 = deltaIdeal*deltaIdeal; float deltaIdeal2 = deltaIdeal*deltaIdeal;
energy += bond.y * deltaIdeal2*( 1.0f + bond.z*deltaIdeal + bond.w*deltaIdeal2 ); energy += bond.y * deltaIdeal2*( 1.0f + cAmoebaSim.amoebaBondCubicParameter*deltaIdeal +
float dEdR = 2.0f*bond.y * deltaIdeal*( 1.0f + 1.5f*bond.z*deltaIdeal + 2.0f*bond.w*deltaIdeal2 ); cAmoebaSim.amoebaBondQuarticicParameter*deltaIdeal2 );
float dEdR = 2.0f*bond.y * deltaIdeal*( 1.0f + 1.5f*cAmoebaSim.amoebaBondCubicParameter*deltaIdeal +
2.0f*cAmoebaSim.amoebaBondQuarticicParameter*deltaIdeal2 );
dEdR = (r > 0.0f) ? (dEdR / r) : 0.0f; dEdR = (r > 0.0f) ? (dEdR / r) : 0.0f;
#endif #endif
......
...@@ -784,9 +784,9 @@ static int readAmoebaHarmonicBondParameters( FILE* filePtr, MapStringInt& forceM ...@@ -784,9 +784,9 @@ static int readAmoebaHarmonicBondParameters( FILE* filePtr, MapStringInt& forceM
methodName.c_str(), arraySize, bondForce->getAmoebaGlobalHarmonicBondCubic(), bondForce->getAmoebaGlobalHarmonicBondQuartic() ); methodName.c_str(), arraySize, bondForce->getAmoebaGlobalHarmonicBondCubic(), bondForce->getAmoebaGlobalHarmonicBondQuartic() );
for( unsigned int ii = 0; ii < arraySize; ii++ ){ for( unsigned int ii = 0; ii < arraySize; ii++ ){
int particle1, particle2; int particle1, particle2;
double length, k, cubic, quartic; double length, k;
bondForce->getBondParameters( ii, particle1, particle2, length, k, cubic, quartic ); bondForce->getBondParameters( ii, particle1, particle2, length, k );
(void) fprintf( log, "%8d %8d %8d %14.7e %14.7e %14.7e %14.7e\n", ii, particle1, particle2, length, k, cubic, quartic); (void) fprintf( log, "%8d %8d %8d %14.7e %14.7e\n", ii, particle1, particle2, length, k );
// skip to end // skip to end
...@@ -799,15 +799,18 @@ static int readAmoebaHarmonicBondParameters( FILE* filePtr, MapStringInt& forceM ...@@ -799,15 +799,18 @@ static int readAmoebaHarmonicBondParameters( FILE* filePtr, MapStringInt& forceM
// convert units to kJ-nm from kCal-Angstrom? // convert units to kJ-nm from kCal-Angstrom?
if( useOpenMMUnits ){ if( useOpenMMUnits ){
double cubic = bondForce->getAmoebaGlobalHarmonicBondCubic()/AngstromToNm;
double quartic = bondForce->getAmoebaGlobalHarmonicBondQuartic()/AngstromToNm*AngstromToNm;
bondForce->setAmoebaGlobalHarmonicBondCubic( cubic );
bondForce->setAmoebaGlobalHarmonicBondQuartic( quartic );
for( int ii = 0; ii < bondForce->getNumBonds(); ii++ ){ for( int ii = 0; ii < bondForce->getNumBonds(); ii++ ){
int particle1, particle2; int particle1, particle2;
double length, k, cubic, quartic; double length, k;
bondForce->getBondParameters( ii, particle1, particle2, length, k, cubic, quartic ); bondForce->getBondParameters( ii, particle1, particle2, length, k );
length *= AngstromToNm; length *= AngstromToNm;
k *= CalToJoule/(AngstromToNm*AngstromToNm); k *= CalToJoule/(AngstromToNm*AngstromToNm);
cubic /= AngstromToNm; bondForce->setBondParameters( ii, particle1, particle2, length, k );
quartic /= AngstromToNm*AngstromToNm;
bondForce->setBondParameters( ii, particle1, particle2, length, k, cubic, quartic );
} }
if( log ){ if( log ){
static const unsigned int maxPrint = MAX_PRINT; static const unsigned int maxPrint = MAX_PRINT;
...@@ -816,8 +819,8 @@ static int readAmoebaHarmonicBondParameters( FILE* filePtr, MapStringInt& forceM ...@@ -816,8 +819,8 @@ static int readAmoebaHarmonicBondParameters( FILE* filePtr, MapStringInt& forceM
methodName.c_str(), arraySize, bondForce->getAmoebaGlobalHarmonicBondCubic(), bondForce->getAmoebaGlobalHarmonicBondQuartic() ); methodName.c_str(), arraySize, bondForce->getAmoebaGlobalHarmonicBondCubic(), bondForce->getAmoebaGlobalHarmonicBondQuartic() );
for( unsigned int ii = 0; ii < arraySize; ii++ ){ for( unsigned int ii = 0; ii < arraySize; ii++ ){
int particle1, particle2; int particle1, particle2;
double length, k, cubic, quartic; double length, k;
bondForce->getBondParameters( ii, particle1, particle2, length, k, cubic, quartic ); bondForce->getBondParameters( ii, particle1, particle2, length, k );
(void) fprintf( log, "%8d %8d %8d %14.7e %14.7e %14.7e %14.7e\n", ii, particle1, particle2, length, k, cubic, quartic); (void) fprintf( log, "%8d %8d %8d %14.7e %14.7e %14.7e %14.7e\n", ii, particle1, particle2, length, k, cubic, quartic);
// skip to end // skip to end
...@@ -944,7 +947,7 @@ static int readAmoebaHarmonicAngleParameters( FILE* filePtr, MapStringInt& force ...@@ -944,7 +947,7 @@ static int readAmoebaHarmonicAngleParameters( FILE* filePtr, MapStringInt& force
for( unsigned int ii = 0; ii < arraySize; ii++ ){ for( unsigned int ii = 0; ii < arraySize; ii++ ){
int particle1, particle2, particle3; int particle1, particle2, particle3;
double length, k, cubic, quartic, pentic, sextic; double length, k;
angleForce->getAngleParameters( ii, particle1, particle2, particle3, length, k ); angleForce->getAngleParameters( ii, particle1, particle2, particle3, length, k );
(void) fprintf( log, "%8d %8d %8d %8d %14.7e %14.7e\n", (void) fprintf( log, "%8d %8d %8d %8d %14.7e %14.7e\n",
ii, particle1, particle2, particle3, length, k ); ii, particle1, particle2, particle3, length, k );
...@@ -962,7 +965,7 @@ static int readAmoebaHarmonicAngleParameters( FILE* filePtr, MapStringInt& force ...@@ -962,7 +965,7 @@ static int readAmoebaHarmonicAngleParameters( FILE* filePtr, MapStringInt& force
if( useOpenMMUnits ){ if( useOpenMMUnits ){
for( int ii = 0; ii < angleForce->getNumAngles(); ii++ ){ for( int ii = 0; ii < angleForce->getNumAngles(); ii++ ){
int particle1, particle2, particle3; int particle1, particle2, particle3;
double length, k, cubic, quartic, pentic, sextic; double length, k;
angleForce->getAngleParameters( ii, particle1, particle2, particle3, length, k ); angleForce->getAngleParameters( ii, particle1, particle2, particle3, length, k );
k *= CalToJoule; k *= CalToJoule;
angleForce->setAngleParameters( ii, particle1, particle2, particle3, length, k ); angleForce->setAngleParameters( ii, particle1, particle2, particle3, length, k );
...@@ -977,7 +980,7 @@ static int readAmoebaHarmonicAngleParameters( FILE* filePtr, MapStringInt& force ...@@ -977,7 +980,7 @@ static int readAmoebaHarmonicAngleParameters( FILE* filePtr, MapStringInt& force
for( unsigned int ii = 0; ii < arraySize; ii++ ){ for( unsigned int ii = 0; ii < arraySize; ii++ ){
int particle1, particle2, particle3; int particle1, particle2, particle3;
double length, k, cubic, quartic, pentic, sextic; double length, k;
angleForce->getAngleParameters( ii, particle1, particle2, particle3, length, k ); angleForce->getAngleParameters( ii, particle1, particle2, particle3, length, k );
(void) fprintf( log, "%8d %8d %8d %8d %14.7e %14.7e\n", (void) fprintf( log, "%8d %8d %8d %8d %14.7e %14.7e\n",
ii, particle1, particle2, particle3, length, k ); ii, particle1, particle2, particle3, length, k );
...@@ -1269,7 +1272,7 @@ static int readAmoebaTorsionParameters( FILE* filePtr, MapStringInt& forceMap, c ...@@ -1269,7 +1272,7 @@ static int readAmoebaTorsionParameters( FILE* filePtr, MapStringInt& forceMap, c
// convert units to kJ-nm from kCal-Angstrom? // convert units to kJ-nm from kCal-Angstrom?
if( useOpenMMUnits ){ if( useOpenMMUnits ){
for( unsigned int ii = 0; ii < torsionForce->getNumTorsions(); ii++ ){ for( int ii = 0; ii < torsionForce->getNumTorsions(); ii++ ){
int particle1, particle2, particle3, particle4; int particle1, particle2, particle3, particle4;
std::vector<double> torsion1; std::vector<double> torsion1;
std::vector<double> torsion2; std::vector<double> torsion2;
...@@ -1405,7 +1408,7 @@ static int readAmoebaPiTorsionParameters( FILE* filePtr, MapStringInt& forceMap, ...@@ -1405,7 +1408,7 @@ static int readAmoebaPiTorsionParameters( FILE* filePtr, MapStringInt& forceMap,
} }
if( useOpenMMUnits ){ if( useOpenMMUnits ){
for( unsigned int ii = 0; ii < piTorsionForce->getNumPiTorsions(); ii++ ){ for( int ii = 0; ii < piTorsionForce->getNumPiTorsions(); ii++ ){
int particle1, particle2, particle3, particle4, particle5, particle6; int particle1, particle2, particle3, particle4, particle5, particle6;
double torsionK; double torsionK;
piTorsionForce->getPiTorsionParameters( ii, particle1, particle2, particle3, particle4, particle5, particle6, torsionK ); piTorsionForce->getPiTorsionParameters( ii, particle1, particle2, particle3, particle4, particle5, particle6, torsionK );
...@@ -1535,7 +1538,7 @@ static int readAmoebaStretchBendParameters( FILE* filePtr, MapStringInt& forceMa ...@@ -1535,7 +1538,7 @@ static int readAmoebaStretchBendParameters( FILE* filePtr, MapStringInt& forceMa
} }
if( useOpenMMUnits ){ if( useOpenMMUnits ){
for( unsigned int ii = 0; ii < stretchBendForce->getNumStretchBends(); ii++ ){ for( int ii = 0; ii < stretchBendForce->getNumStretchBends(); ii++ ){
int particle1, particle2, particle3; int particle1, particle2, particle3;
double lengthAB, lengthCB, angle, k; double lengthAB, lengthCB, angle, k;
stretchBendForce->getStretchBendParameters( ii, particle1, particle2, particle3, lengthAB, lengthCB, angle, k ); stretchBendForce->getStretchBendParameters( ii, particle1, particle2, particle3, lengthAB, lengthCB, angle, k );
...@@ -1704,7 +1707,7 @@ static int readAmoebaOutOfPlaneBendParameters( FILE* filePtr, MapStringInt& forc ...@@ -1704,7 +1707,7 @@ static int readAmoebaOutOfPlaneBendParameters( FILE* filePtr, MapStringInt& forc
} }
if( useOpenMMUnits ){ if( useOpenMMUnits ){
for( unsigned int ii = 0; ii < outOfPlaneBendForce->getNumOutOfPlaneBends(); ii++ ){ for( int ii = 0; ii < outOfPlaneBendForce->getNumOutOfPlaneBends(); ii++ ){
int particle1, particle2, particle3, particle4; int particle1, particle2, particle3, particle4;
double k; double k;
outOfPlaneBendForce->getOutOfPlaneBendParameters( ii, particle1, particle2, particle3, particle4, k ); outOfPlaneBendForce->getOutOfPlaneBendParameters( ii, particle1, particle2, particle3, particle4, k );
...@@ -1776,7 +1779,7 @@ static int readAmoebaTorsionTorsionGrid( FILE* filePtr, int numX, int numY, Tors ...@@ -1776,7 +1779,7 @@ static int readAmoebaTorsionTorsionGrid( FILE* filePtr, int numX, int numY, Tors
} }
grid.resize( numX ); grid.resize( numX );
for( unsigned int ii = 0; ii < numX; ii++ ){ for( int ii = 0; ii < numX; ii++ ){
grid[ii].resize( numY ); grid[ii].resize( numY );
} }
for( int ii = 0; ii < gridCount; ii++ ){ for( int ii = 0; ii < gridCount; ii++ ){
...@@ -1821,7 +1824,7 @@ static int readAmoebaTorsionTorsionGrid( FILE* filePtr, int numX, int numY, Tors ...@@ -1821,7 +1824,7 @@ static int readAmoebaTorsionTorsionGrid( FILE* filePtr, int numX, int numY, Tors
// 'top' of grid // 'top' of grid
for( unsigned int ii = 0; ii < gridCount; ii++ ){ for( int ii = 0; ii < gridCount; ii++ ){
std::vector<double> values = grid[xI][yI]; std::vector<double> values = grid[xI][yI];
(void) fprintf( log, "%4d %4d %4d ", ii, xI, yI ); (void) fprintf( log, "%4d %4d %4d ", ii, xI, yI );
for( unsigned int jj = 0; jj < values.size(); jj++ ){ for( unsigned int jj = 0; jj < values.size(); jj++ ){
...@@ -1843,7 +1846,7 @@ static int readAmoebaTorsionTorsionGrid( FILE* filePtr, int numX, int numY, Tors ...@@ -1843,7 +1846,7 @@ static int readAmoebaTorsionTorsionGrid( FILE* filePtr, int numX, int numY, Tors
xI = 0; xI = 0;
yI = (gridCount/numX) - 3; yI = (gridCount/numX) - 3;
for( unsigned int ii = 0; ii < gridCount; ii++ ){ for( int ii = 0; ii < gridCount; ii++ ){
std::vector<double> values = grid[xI][yI]; std::vector<double> values = grid[xI][yI];
(void) fprintf( log, "%4d %4d %4d ", (void) fprintf( log, "%4d %4d %4d ",
ii, xI, yI ); ii, xI, yI );
...@@ -2038,7 +2041,7 @@ static void readAmoebaMultipoleCovalent( FILE* filePtr, AmoebaMultipoleForce* mu ...@@ -2038,7 +2041,7 @@ static void readAmoebaMultipoleCovalent( FILE* filePtr, AmoebaMultipoleForce* mu
char* isNotEof = "1"; char* isNotEof = "1";
for( int ii = 0; ii < numberOfMultipoles && isNotEof; ii++ ){ for( int ii = 0; ii < numberOfMultipoles && isNotEof; ii++ ){
for( int jj = 0; jj < numberOfCovalentTypes && isNotEof; jj++ ){ for( unsigned int jj = 0; jj < numberOfCovalentTypes && isNotEof; jj++ ){
StringVector lineTokens; StringVector lineTokens;
isNotEof = readLine( filePtr, lineTokens, lineCount, log ); isNotEof = readLine( filePtr, lineTokens, lineCount, log );
std::vector<int> intVector; std::vector<int> intVector;
...@@ -2273,12 +2276,12 @@ scalingDistanceCutoff, thole, dampingFactor, pGamma, ...@@ -2273,12 +2276,12 @@ scalingDistanceCutoff, thole, dampingFactor, pGamma,
double polarityConversion = AngstromToNm*AngstromToNm*AngstromToNm; double polarityConversion = AngstromToNm*AngstromToNm*AngstromToNm;
double dampingFactorConversion = sqrt( AngstromToNm ); double dampingFactorConversion = sqrt( AngstromToNm );
float scalingDistanceCutoff = multipoleForce->getScalingDistanceCutoff(); float scalingDistanceCutoff = static_cast<float>(multipoleForce->getScalingDistanceCutoff());
scalingDistanceCutoff *= AngstromToNm; scalingDistanceCutoff *= static_cast<float>(AngstromToNm);
multipoleForce->setScalingDistanceCutoff( scalingDistanceCutoff ); multipoleForce->setScalingDistanceCutoff( scalingDistanceCutoff );
for( unsigned int ii = 0; ii < multipoleForce->getNumMultipoles(); ii++ ){ for( int ii = 0; ii < multipoleForce->getNumMultipoles(); ii++ ){
int axisType, zAxis, xAxis; int axisType, zAxis, xAxis;
std::vector<double> dipole; std::vector<double> dipole;
...@@ -2299,8 +2302,8 @@ scalingDistanceCutoff, thole, dampingFactor, pGamma, ...@@ -2299,8 +2302,8 @@ scalingDistanceCutoff, thole, dampingFactor, pGamma,
} }
} else { } else {
float electricConstant = multipoleForce->getElectricConstant(); float electricConstant = static_cast<float>(multipoleForce->getElectricConstant());
electricConstant /= (AngstromToNm*CalToJoule); electricConstant /= static_cast<float>(AngstromToNm*CalToJoule);
multipoleForce->setElectricConstant( electricConstant ); multipoleForce->setElectricConstant( electricConstant );
} }
...@@ -2475,7 +2478,7 @@ static int readAmoebaGeneralizedKirkwoodParameters( FILE* filePtr, MapStringInt& ...@@ -2475,7 +2478,7 @@ static int readAmoebaGeneralizedKirkwoodParameters( FILE* filePtr, MapStringInt&
} }
if( useOpenMMUnits ){ if( useOpenMMUnits ){
for( unsigned int ii = 0; ii < gbsaObcForce->getNumParticles(); ii++ ){ for( int ii = 0; ii < gbsaObcForce->getNumParticles(); ii++ ){
double charge, radius, scalingFactor; double charge, radius, scalingFactor;
gbsaObcForce->getParticleParameters( ii, charge, radius, scalingFactor ); gbsaObcForce->getParticleParameters( ii, charge, radius, scalingFactor );
radius *= AngstromToNm; radius *= AngstromToNm;
...@@ -2683,7 +2686,7 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const ...@@ -2683,7 +2686,7 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const
if( log ){ if( log ){
//static const unsigned int maxPrint = MAX_PRINT; //static const unsigned int maxPrint = MAX_PRINT;
static const unsigned int maxPrint = 15; static const int maxPrint = 15;
unsigned int arraySize = static_cast<unsigned int>(vdwForce->getNumParticles()); unsigned int arraySize = static_cast<unsigned int>(vdwForce->getNumParticles());
unsigned int tableSize = static_cast<unsigned int>(vdwForce->getSigEpsTableSize()); unsigned int tableSize = static_cast<unsigned int>(vdwForce->getSigEpsTableSize());
(void) fprintf( log, "%s: %u sample of AmoebaVdwForce parameters; SigEpsTableSize=%d combining rules=[sig=%s eps=%s]\n", (void) fprintf( log, "%s: %u sample of AmoebaVdwForce parameters; SigEpsTableSize=%d combining rules=[sig=%s eps=%s]\n",
...@@ -2714,7 +2717,7 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const ...@@ -2714,7 +2717,7 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const
(void) fprintf( log, "SigWEps table not loaded\n" ); (void) fprintf( log, "SigWEps table not loaded\n" );
} }
for( unsigned int ii = 0; ii < vdwForce->getNumParticles(); ii++ ){ for( int ii = 0; ii < vdwForce->getNumParticles(); ii++ ){
int indexIV, indexClass; int indexIV, indexClass;
double sigma, sigma4, epsilon, epsilon4, reduction; double sigma, sigma4, epsilon, epsilon4, reduction;
std::vector< int > exclusions; std::vector< int > exclusions;
...@@ -2731,7 +2734,7 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const ...@@ -2731,7 +2734,7 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const
// skip to end // skip to end
if( ii == maxPrint && (arraySize - maxPrint) > ii ){ if( ii == maxPrint && static_cast<int>(arraySize - maxPrint) > ii ){
ii = arraySize - maxPrint - 1; ii = arraySize - maxPrint - 1;
} }
} }
...@@ -2791,7 +2794,7 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const ...@@ -2791,7 +2794,7 @@ static int readAmoebaVdwParameters( FILE* filePtr, MapStringInt& forceMap, const
*/ */
} }
for( unsigned int ii = 0; ii < vdwForce->getNumParticles(); ii++ ){ for( int ii = 0; ii < vdwForce->getNumParticles(); ii++ ){
int indexIV, indexClass; int indexIV, indexClass;
double sigma, sigma4, epsilon, epsilon4, reduction; double sigma, sigma4, epsilon, epsilon4, reduction;
vdwForce->getParticleParameters( ii, indexIV, indexClass, sigma, sigma4, epsilon, epsilon4, reduction ); vdwForce->getParticleParameters( ii, indexIV, indexClass, sigma, sigma4, epsilon, epsilon4, reduction );
...@@ -2873,7 +2876,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force ...@@ -2873,7 +2876,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force
double epsilon = atof( lineTokens[tokenIndex++].c_str() ); double epsilon = atof( lineTokens[tokenIndex++].c_str() );
wcaDispersionForce->addParticle( radius, epsilon ); wcaDispersionForce->addParticle( radius, epsilon );
if( tokenIndex < lineTokens.size() ){ if( tokenIndex < static_cast<int>(lineTokens.size()) ){
double cdisp = atof( lineTokens[tokenIndex++].c_str() ); double cdisp = atof( lineTokens[tokenIndex++].c_str() );
maxDispersionEnergyVector.push_back( cdisp ); maxDispersionEnergyVector.push_back( cdisp );
} }
...@@ -3015,7 +3018,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force ...@@ -3015,7 +3018,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force
wcaDispersionForce->setEpso( epso ); wcaDispersionForce->setEpso( epso );
wcaDispersionForce->setEpsh( epsh ); wcaDispersionForce->setEpsh( epsh );
for( unsigned int ii = 0; ii < wcaDispersionForce->getNumParticles(); ii++ ){ for( int ii = 0; ii < wcaDispersionForce->getNumParticles(); ii++ ){
double radius, epsilon, maxDispersionEnergy; double radius, epsilon, maxDispersionEnergy;
wcaDispersionForce->getParticleParameters( ii, radius, epsilon ); wcaDispersionForce->getParticleParameters( ii, radius, epsilon );
...@@ -3023,7 +3026,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force ...@@ -3023,7 +3026,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force
epsilon *= CalToJoule; epsilon *= CalToJoule;
wcaDispersionForce->setParticleParameters( ii, radius, epsilon ); wcaDispersionForce->setParticleParameters( ii, radius, epsilon );
if( ii < maxDispersionEnergyVector.size() ){ if( ii < static_cast<int>(maxDispersionEnergyVector.size()) ){
wcaDispersionForce->getMaximumDispersionEnergy( ii, maxDispersionEnergy ); wcaDispersionForce->getMaximumDispersionEnergy( ii, maxDispersionEnergy );
double tinkerValue = maxDispersionEnergyVector[ii]; double tinkerValue = maxDispersionEnergyVector[ii];
tinkerValue *= CalToJoule; tinkerValue *= CalToJoule;
...@@ -3694,7 +3697,6 @@ Integrator* readAmoebaParameterFile( const std::string& inputParameterFile, MapS ...@@ -3694,7 +3697,6 @@ Integrator* readAmoebaParameterFile( const std::string& inputParameterFile, MapS
static const std::string methodName = "readAmoebaParameterFile"; static const std::string methodName = "readAmoebaParameterFile";
int PrintOn = 1; int PrintOn = 1;
double unitsConversion[LastUnits];
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -4013,7 +4015,7 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe ...@@ -4013,7 +4015,7 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe
int useOpenMMUnits, FILE* log ) { int useOpenMMUnits, FILE* log ) {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
/*
static const std::string methodName = "checkIntermediateMultipoleQuantities"; static const std::string methodName = "checkIntermediateMultipoleQuantities";
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -4074,7 +4076,7 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe ...@@ -4074,7 +4076,7 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe
totalMisses += misses; totalMisses += misses;
} }
} catch( exception& e ){ } catch( exception& e ){
(void) fprintf( log, "Rotation matricies not available\n" ); (void) fprintf( log, "Rotation matricies not available %s\n", e.what() );
(void) fflush( log ); (void) fflush( log );
} }
...@@ -4122,7 +4124,7 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe ...@@ -4122,7 +4124,7 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe
totalMisses += misses; totalMisses += misses;
} }
} catch( exception& e ){ } catch( exception& e ){
(void) fprintf( log, "Fixed-E fields not available\n" ); (void) fprintf( log, "Fixed-E fields not available %s\n", e.what() );
(void) fflush( log ); (void) fflush( log );
} }
...@@ -4172,10 +4174,10 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe ...@@ -4172,10 +4174,10 @@ void checkIntermediateMultipoleQuantities( Context* context, MapStringVectorOfVe
totalMisses += misses; totalMisses += misses;
} }
} catch( exception& e ){ } catch( exception& e ){
(void) fprintf( log, "Induced dipoles not available\n" ); (void) fprintf( log, "Induced dipoles not available %s\n", e.what() );
(void) fflush( log ); (void) fflush( log );
} }
*/
} }
void calculateBorn1( System& amoebaSystem, std::vector<Vec3>& tinkerCoordinates, FILE* log ) { void calculateBorn1( System& amoebaSystem, std::vector<Vec3>& tinkerCoordinates, FILE* log ) {
...@@ -4636,7 +4638,7 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete ...@@ -4636,7 +4638,7 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete
std::vector<Vec3> expectedForces; std::vector<Vec3> expectedForces;
expectedForces.resize( system.getNumParticles() ); expectedForces.resize( system.getNumParticles() );
for( unsigned int ii = 0; ii < system.getNumParticles(); ii++ ){ for( int ii = 0; ii < system.getNumParticles(); ii++ ){
expectedForces[ii][0] = 0.0; expectedForces[ii][0] = 0.0;
expectedForces[ii][1] = 0.0; expectedForces[ii][1] = 0.0;
expectedForces[ii][2] = 0.0; expectedForces[ii][2] = 0.0;
...@@ -4646,7 +4648,7 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete ...@@ -4646,7 +4648,7 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete
for( unsigned int ii = 0; ii < forceList.size(); ii++ ){ for( unsigned int ii = 0; ii < forceList.size(); ii++ ){
expectedEnergy += tinkerEnergies[forceList[ii]]; expectedEnergy += tinkerEnergies[forceList[ii]];
std::vector<Vec3> forces = tinkerForces[forceList[ii]]; std::vector<Vec3> forces = tinkerForces[forceList[ii]];
for( unsigned int jj = 0; jj < system.getNumParticles(); jj++ ){ for( int jj = 0; jj < system.getNumParticles(); jj++ ){
expectedForces[jj][0] += forces[jj][0]; expectedForces[jj][0] += forces[jj][0];
expectedForces[jj][1] += forces[jj][1]; expectedForces[jj][1] += forces[jj][1];
expectedForces[jj][2] += forces[jj][2]; expectedForces[jj][2] += forces[jj][2];
...@@ -4757,6 +4759,10 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete ...@@ -4757,6 +4759,10 @@ void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParamete
} }
} }
int isNanOrInfinity( double number ){
return (number != number || number == std::numeric_limits<double>::infinity() || number == -std::numeric_limits<double>::infinity()) ? 1 : 0;
}
/** /**
* Check that energy and force are consistent * Check that energy and force are consistent
* *
...@@ -4825,15 +4831,15 @@ void testEnergyForcesConsistent( std::string parameterFileName, MapStringInt& fo ...@@ -4825,15 +4831,15 @@ void testEnergyForcesConsistent( std::string parameterFileName, MapStringInt& fo
// check norm is not nan // check norm is not nan
if( isinf( forceNorm ) || isnan( forceNorm ) ){ if( isNanOrInfinity( forceNorm ) ){
if( log ){ if( log ){
(void) fprintf( log, "%s norm of force is nan -- aborting.\n", methodName.c_str() ); (void) fprintf( log, "%s norm of force is nan -- aborting.\n", methodName.c_str() );
unsigned int hitNan = 0; unsigned int hitNan = 0;
for( unsigned int ii = 0; (ii < forces.size()) && (hitNan < 10); ii++ ){ for( unsigned int ii = 0; (ii < forces.size()) && (hitNan < 10); ii++ ){
if( isinf( forces[ii][0] ) || isnan( forces[ii][0] ) || if( isNanOrInfinity( forces[ii][0] ) ||
isinf( forces[ii][1] ) || isnan( forces[ii][1] ) || isNanOrInfinity( forces[ii][1] ) ||
isinf( forces[ii][2] ) || isnan( forces[ii][2] ) )hitNan++; isNanOrInfinity( forces[ii][2] ) )hitNan++;
(void) fprintf( log, "%6u x[%14.7e %14.7e %14.7e] f[%14.7e %14.7e %14.7e]\n", ii, (void) fprintf( log, "%6u x[%14.7e %14.7e %14.7e] f[%14.7e %14.7e %14.7e]\n", ii,
coordinates[ii][0], coordinates[ii][1], coordinates[ii][2], coordinates[ii][0], coordinates[ii][1], coordinates[ii][2],
...@@ -5121,7 +5127,7 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM ...@@ -5121,7 +5127,7 @@ void testEnergyConservation( std::string parameterFileName, MapStringInt& forceM
FILE* intermediateStateFile = NULL; FILE* intermediateStateFile = NULL;
if( intermediateStateFileName != "NA" ){ if( intermediateStateFileName != "NA" ){
#ifdef _MSC_VER #ifdef _MSC_VER
fopen_s( &intermediateStateFile, intermediateStateFileName.c_str(), "w")); fopen_s( &intermediateStateFile, intermediateStateFileName.c_str(), "w");
#else #else
intermediateStateFile = fopen( intermediateStateFileName.c_str(), "w" ); intermediateStateFile = fopen( intermediateStateFileName.c_str(), "w" );
#endif #endif
......
...@@ -48,6 +48,7 @@ ...@@ -48,6 +48,7 @@
#include "AmoebaVdwForce.h" #include "AmoebaVdwForce.h"
#include "AmoebaWcaDispersionForce.h" #include "AmoebaWcaDispersionForce.h"
#include "AmoebaSASAForce.h" #include "AmoebaSASAForce.h"
#include "internal/windowsExport.h"
#include <ctime> #include <ctime>
#include <vector> #include <vector>
...@@ -232,5 +233,6 @@ void initializeForceMap( MapStringInt& forceMap, int initialValue ); ...@@ -232,5 +233,6 @@ void initializeForceMap( MapStringInt& forceMap, int initialValue );
void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParameterFileName, MapStringInt& forceMap, void testUsingAmoebaTinkerParameterFile( const std::string& amoebaTinkerParameterFileName, MapStringInt& forceMap,
double tolerance, FILE* summaryFile, FILE* log ); double tolerance, FILE* summaryFile, FILE* log );
int runTestsUsingAmoebaTinkerParameterFile( MapStringString& argumentMap );
void appendInputArgumentsToArgumentMap( int numberOfArguments, char* argv[], MapStringString& argumentMap ); int OPENMM_EXPORT runTestsUsingAmoebaTinkerParameterFile( MapStringString& argumentMap );
void OPENMM_EXPORT appendInputArgumentsToArgumentMap( int numberOfArguments, char* argv[], MapStringString& argumentMap );
...@@ -5,20 +5,26 @@ ENABLE_TESTING() ...@@ -5,20 +5,26 @@ ENABLE_TESTING()
# INCLUDE(${CMAKE_SOURCE_DIR}/platforms/cuda/cuda-cmake/FindCuda.cmake) # INCLUDE(${CMAKE_SOURCE_DIR}/platforms/cuda/cuda-cmake/FindCuda.cmake)
INCLUDE_DIRECTORIES(${CUDA_INCLUDE}) INCLUDE_DIRECTORIES(${CUDA_INCLUDE})
INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/include) INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/include)
INCLUDE_DIRECTORIES(${OPENMM_DIR}/openmmapi/include/openmm)
INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/src) INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/src)
INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/src/kernels) INCLUDE_DIRECTORIES(${OPENMM_DIR}/platforms/cuda/src/kernels)
SET(SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET "AmoebaTinkerParameterFile" ) SET(SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET "AmoebaTinkerParameterFile" )
SET(AMOEBA_TINKER_PARAMETER_FILE_SOURCE_FILES "AmoebaTinkerParameterFile.cpp" ) SET(AMOEBA_TINKER_PARAMETER_FILE_SOURCE_FILES "AmoebaTinkerParameterFile.cpp" )
SET(AMOEBA_TINKER_PARAMETER_FILE_INCLUDE_FILES "AmoebaTinkerParameterFile.h" ) SET(AMOEBA_TINKER_PARAMETER_FILE_INCLUDE_FILES "AmoebaTinkerParameterFile.h" )
ADD_LIBRARY(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} SHARED ${AMOEBA_TINKER_PARAMETER_FILE_SOURCE_FILES} ${AMOEBA_TINKER_PARAMETER_FILE_INCLUDE_FILES} )
SET_TARGET_PROPERTIES(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_BUILDING_SHARED_LIBRARY -DLEPTON_BUILDING_SHARED_LIBRARY -DOPENMM_VALIDATE_BUILDING_SHARED_LIBRARY")
Set( SHARED_OPENMM__AMOEBA_TARGET OpenMMAmoeba) Set( SHARED_OPENMM__AMOEBA_TARGET OpenMMAmoeba)
Set( STATIC_OPENMM_TARGET OpenMMAmoeba_static) Set( STATIC_OPENMM_TARGET OpenMMAmoeba_static)
Set( SHARED_CUDA_TARGET OpenMMAmoebaCuda ) Set( SHARED_CUDA_TARGET OpenMMAmoebaCuda )
Set( STATIC_CUDA_TARGET OpenMMCuda_static OpenMMAmoebaCuda_static) Set( STATIC_CUDA_TARGET OpenMMCuda_static OpenMMAmoebaCuda_static)
ADD_LIBRARY(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} SHARED ${AMOEBA_TINKER_PARAMETER_FILE_SOURCE_FILES} ${AMOEBA_TINKER_PARAMETER_FILE_INCLUDE_FILES} )
SET_TARGET_PROPERTIES(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_BUILDING_SHARED_LIBRARY -DLEPTON_BUILDING_SHARED_LIBRARY -DOPENMM_VALIDATE_BUILDING_SHARED_LIBRARY")
TARGET_LINK_LIBRARIES(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} debug ${OPENMM_LIBRARY_NAME}_d optimized ${OPENMM_LIBRARY_NAME})
TARGET_LINK_LIBRARIES(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} debug ${SHARED_OPENMM__AMOEBA_TARGET}_d optimized ${SHARED_OPENMM__AMOEBA_TARGET})
TARGET_LINK_LIBRARIES(${SHARED_AMOEBA_TINKER_PARAMETER_FILE_TARGET} debug ${SHARED_CUDA_TARGET}_d optimized ${SHARED_CUDA_TARGET})
IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug) IF (UNIX AND CMAKE_BUILD_TYPE MATCHES Debug)
SET(SHARED_CUDA_TARGET ${SHARED_CUDA_TARGET}_d) SET(SHARED_CUDA_TARGET ${SHARED_CUDA_TARGET}_d)
SET(SHARED_OPENMM__AMOEBA_TARGET ${SHARED_OPENMM__AMOEBA_TARGET}_d) SET(SHARED_OPENMM__AMOEBA_TARGET ${SHARED_OPENMM__AMOEBA_TARGET}_d)
......
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