Skip to content

GitLab

Commit 0cae89d3 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Mods for wrapper code; added yAxis parameter to multipole force-field

parent 077a93c8
Hide whitespace changes
Inline Side-by-side
Showing with 300 additions and 251 deletions
plugins/amoeba/CMakeLists.txt
View file @ 0cae89d3
......@@ -16,24 +16,6 @@
#INCLUDE(Dart)
SET(LOG FALSE)
# ----------------------------------------------------------------------------
IF(LOG)
SET(LOG_FILE "CMakeLog.txt" )
FILE( WRITE ${LOG_FILE} "In OpenMM Cmake\n")
ENDIF(LOG)
IF(LOG)
MACRO(LOG_DIR LOG_FILE DIR_LIST )
FILE( APPEND ${LOG_FILE} "\n${DIR_LIST}\n")
FOREACH(currentFile ${ARGN})
FILE( APPEND ${LOG_FILE} " ${currentFile}\n" )
ENDFOREACH(currentFile)
ENDMACRO(LOG_DIR)
ENDIF(LOG)
# ----------------------------------------------------------------------------
# The source is organized into subdirectories, but we handle them all from
......@@ -107,23 +89,15 @@ INCLUDE_DIRECTORIES(BEFORE ${OPENMM_DIR}/platforms/reference/src/SimTKReference)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/platforms/reference/src)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/platforms/reference/src/SimTKReference)
# ----------------------------------------------------------------------------
IF(LOG)
FILE( APPEND ${LOG_FILE} "SOURCE_AMOEBA_FILES=" ${SOURCE_AMOEBA_FILES} "\n" )
FILE( APPEND ${LOG_FILE} "SOURCE_AMOEBA_INCLUDE_FILES=" ${SOURCE_AMOEBA_INCLUDE_FILES} "\n" )
FILE( APPEND ${LOG_FILE} "OPENMM_PATH=" ${OPENMM_PATH} "\n" )
FILE( APPEND ${LOG_FILE} "CUDA_PATH=" ${CUDA_PATH} "\n" )
ENDIF(LOG)
# ----------------------------------------------------------------------------
# If API_AMOEBA wrappers are being generated, and add them to the build.
#SET(OPENMM_BUILD_API_AMOEBA_WRAPPERS OFF CACHE BOOL "Build wrappers for C and Fortran")
#IF(OPENMM_BUILD_API_AMOEBA_WRAPPERS)
# ADD_SUBDIRECTORY(wrappers)
# SET(SOURCE_AMOEBA_FILES ${SOURCE_AMOEBA_FILES} wrappers/OpenMMCWrapper.cpp wrappers/OpenMMFortranWrapper.cpp)
# SET_SOURCE_AMOEBA_FILES_PROPERTIES(wrappers/OpenMMCWrapper.cpp wrappers/OpenMMFortranWrapper.cpp PROPERTIES GENERATED TRUE)
#ENDIF(OPENMM_BUILD_API_AMOEBA_WRAPPERS)
IF(OPENMM_BUILD_API_WRAPPERS)
ADD_SUBDIRECTORY(wrappers)
SET(SOURCE_AMOEBA_FILES ${SOURCE_AMOEBA_FILES} wrappers/AmoebaOpenMMCWrapper.cpp wrappers/AmoebaOpenMMFortranWrapper.cpp)
SET_SOURCE_FILES_PROPERTIES(wrappers/AmoebaOpenMMCWrapper.cpp wrappers/AmoebaOpenMMFortranWrapper.cpp PROPERTIES GENERATED TRUE)
ENDIF(OPENMM_BUILD_API_WRAPPERS)
INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_SOURCE_DIR}/src)
......@@ -135,10 +109,14 @@ IF(OPENMM_BUILD_STATIC_LIB)
SET_TARGET_PROPERTIES(${STATIC_AMOEBA_TARGET} PROPERTIES COMPILE_FLAGS "-DOPENMM_USE_STATIC_LIBRARIES -DOPENMM_BUILDING_STATIC_LIBRARY -DLEPTON_USE_STATIC_LIBRARIES -DLEPTON_BUILDING_STATIC_LIBRARY")
ENDIF(OPENMM_BUILD_STATIC_LIB)
#IF(OPENMM_BUILD_API_AMOEBA_WRAPPERS)
# ADD_DEPENDENCIES(${SHARED_AMOEBA_TARGET} ApiWrappers)
# ADD_DEPENDENCIES(${STATIC_AMOEBA_TARGET} ApiWrappers)
#ENDIF(OPENMM_BUILD_API_AMOEBA_WRAPPERS)
IF(OPENMM_BUILD_API_WRAPPERS)
ADD_DEPENDENCIES(${SHARED_AMOEBA_TARGET} AmoebaApiWrappers)
IF(OPENMM_BUILD_STATIC_LIB)
ADD_DEPENDENCIES(${STATIC_AMOEBA_TARGET} AmoebaApiWrappers)
ENDIF(OPENMM_BUILD_STATIC_LIB)
ENDIF(OPENMM_BUILD_API_WRAPPERS)
# ----------------------------------------------------------------------------
# On Linux need to link to libdl
FIND_LIBRARY(DL_LIBRARY dl)
......@@ -177,11 +155,6 @@ IF(OPENMM_BUILD_AMOEBA_CUDA_LIB)
SET(OPENMM_AMOEBA_CUDA_SOURCE_SUBDIRS . openmmapi olla platforms/cuda)
ENDIF(OPENMM_BUILD_AMOEBA_CUDA_LIB)
#SET(OPENMM_BUILD_BROOK_LIB OFF CACHE BOOL "Build OpenMMBrook library for ATI GPUs")
#IF(OPENMM_BUILD_BROOK_LIB)
# ADD_SUBDIRECTORY(platforms/brook)
#ENDIF(OPENMM_BUILD_BROOK_LIB)
#
#SET(OPENMM_BUILD_OPENCL_LIB OFF CACHE BOOL "Build OpenMMOpenCL library for Nvidia GPUs")
#IF(OPENMM_BUILD_OPENCL_LIB)
# ADD_SUBDIRECTORY(platforms/opencl)
......@@ -210,7 +183,6 @@ INSTALL_FILES(/include/openmm/internal FILES ${INTERNAL_HEADERS})
# ENDIF (NOT CMAKE_BUILD_TYPE OR CMAKE_BUILD_TYPE MATCHES Debug)
#ENDIF (UNIX AND NOT CYGWIN AND NOT APPLE)
#
# Testing
#
......
plugins/amoeba/openmmapi/include/AmoebaMultipoleForce.h
View file @ 0cae89d3
......@@ -162,8 +162,9 @@ public:
* @param molecularDipole the particle's molecular dipole (vector of size 3)
* @param molecularQuadrupole the particle's molecular quadrupole (vector of size 9)
* @param axisType the particle's axis type ( ZThenX, Bisector )
* @param multipoleAtomId1 index of first atom used in constructing lab<->molecular frames
* @param multipoleAtomId2 index of second atom used in constructing lab<->molecular frames
* @param multipoleAtomZ index of first atom used in constructing lab<->molecular frames
* @param multipoleAtomX index of second atom used in constructing lab<->molecular frames
* @param multipoleAtomY index of second atom used in constructing lab<->molecular frames
* @param thole Thole parameter
* @param dampingFactor dampingFactor parameter
* @param polarity polarity parameter
......@@ -171,7 +172,7 @@ public:
* @return the index of the particle that was added
*/
int addParticle( double charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole, int axisType,
int multipoleAtomId1, int multipoleAtomId2, double thole, double dampingFactor, double polarity );
int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity );
/**
* Get the multipole parameters for a particle.
......@@ -181,14 +182,15 @@ public:
* @param molecularDipole the particle's molecular dipole (vector of size 3)
* @param molecularQuadrupole the particle's molecular quadrupole (vector of size 9)
* @param axisType the particle's axis type ( ZThenX, Bisector )
* @param multipoleAtomId1 index of first atom used in constructing lab<->molecular frames
* @param multipoleAtomId2 index of second atom used in constructing lab<->molecular frames
* @param multipoleAtomZ index of first atom used in constructing lab<->molecular frames
* @param multipoleAtomX index of second atom used in constructing lab<->molecular frames
* @param multipoleAtomY index of second atom used in constructing lab<->molecular frames
* @param thole Thole parameter
* @param dampingFactor dampingFactor parameter
* @param polarity polarity parameter
*/
void getMultipoleParameters(int index, double& charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole,
int& axisType, int& multipoleAtomId1, int& multipoleAtomId2, double& thole, double& dampingFactor, double& polarity ) const;
int& axisType, int& multipoleAtomZ, int& multipoleAtomX, int& multipoleAtomY, double& thole, double& dampingFactor, double& polarity ) const;
/**
* Set the multipole parameters for a particle.
......@@ -198,12 +200,13 @@ public:
* @param molecularDipole the particle's molecular dipole (vector of size 3)
* @param molecularQuadrupole the particle's molecular quadrupole (vector of size 9)
* @param axisType the particle's axis type ( ZThenX, Bisector )
* @param multipoleAtomId1 index of first atom used in constructing lab<->molecular frames
* @param multipoleAtomId2 index of second atom used in constructing lab<->molecular frames
* @param multipoleAtomZ index of first atom used in constructing lab<->molecular frames
* @param multipoleAtomX index of second atom used in constructing lab<->molecular frames
* @param multipoleAtomY index of second atom used in constructing lab<->molecular frames
* @param polarity polarity parameter
*/
void setMultipoleParameters(int index, double charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole,
int axisType, int multipoleAtomId1, int multipoleAtomId2, double thole, double dampingFactor, double polarity);
int axisType, int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity);
/**
* Set the CovalentMap for an atom
......@@ -338,7 +341,7 @@ private:
class AmoebaMultipoleForce::MultipoleInfo {
public:
int axisType, multipoleAtomId1, multipoleAtomId2;
int axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY;
double charge, thole, dampingFactor, polarity;
std::vector<double> molecularDipole;
......@@ -346,7 +349,7 @@ public:
std::vector< std::vector<int> > covalentInfo;
MultipoleInfo() {
axisType = multipoleAtomId1 = multipoleAtomId2 = -1;
axisType = multipoleAtomZ = multipoleAtomX = multipoleAtomY = -1;
charge = thole = dampingFactor = 0.0;
molecularDipole.resize( 3 );
......@@ -355,8 +358,8 @@ public:
}
MultipoleInfo( double charge, std::vector<double>& inputMolecularDipole, std::vector<double>& inputMolecularQuadrupole,
int axisType, int multipoleAtomId1, int multipoleAtomId2, double thole, double dampingFactor, double polarity) :
charge(charge), axisType(axisType), multipoleAtomId1(multipoleAtomId1), multipoleAtomId2(multipoleAtomId2),
int axisType, int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity) :
charge(charge), axisType(axisType), multipoleAtomZ(multipoleAtomZ), multipoleAtomX(multipoleAtomX), multipoleAtomY(multipoleAtomY),
thole(thole), dampingFactor(dampingFactor), polarity(polarity) {
covalentInfo.resize( CovalentEnd );
......
plugins/amoeba/openmmapi/include/AmoebaTorsionForce.h
View file @ 0cae89d3
......@@ -75,7 +75,7 @@ public:
* @return the index of the torsion that was added
*/
int addTorsion(int particle1, int particle2, int particle3, int particle4,
std::vector<double> torsion1, std::vector<double> torsion2, std::vector<double> torsion3 );
std::vector<double>& torsion1, std::vector<double>& torsion2, std::vector<double>& torsion3 );
/**
* Get the force field parameters for a torsion term.
......@@ -105,7 +105,7 @@ public:
* @param torsion3 the vector of torsion params for third index (amplitude, phase, fold)
*/
void setTorsionParameters(int index, int particle1, int particle2, int particle3, int particle4,
std::vector<double> torsion1, std::vector<double> torsion2, std::vector<double> torsion3 );
std::vector<double>& torsion1, std::vector<double>& torsion2, std::vector<double>& torsion3 );
protected:
ForceImpl* createImpl();
......
plugins/amoeba/openmmapi/include/AmoebaTorsionTorsionForce.h
View file @ 0cae89d3
......@@ -117,7 +117,8 @@ public:
* @param gridIndex the grid index
* @return grid return grid reference
*/
const TorsionTorsionGrid& getTorsionTorsionGrid( int index ) const;
//const TorsionTorsionGrid& getTorsionTorsionGrid( int index ) const;
const std::vector< std::vector< std::vector<double> > >& getTorsionTorsionGrid( int index ) const;
/**
* Set the torsion-torsion grid at the specified index
......@@ -131,7 +132,8 @@ public:
* grid[x][y][4] = dfdy value
* grid[x][y][5] = dfd(xy) value
*/
void setTorsionTorsionGrid(int index, TorsionTorsionGrid& grid );
//void setTorsionTorsionGrid(int index, TorsionTorsionGrid& grid );
void setTorsionTorsionGrid(int index, std::vector< std::vector< std::vector<double> > >& grid );
protected:
ForceImpl* createImpl();
......
plugins/amoeba/openmmapi/include/AmoebaVdwForce.h
View file @ 0cae89d3
......@@ -104,28 +104,28 @@ public:
*
* @param sigmaCombiningRule sigma combining rule: 'ARITHMETIC', 'GEOMETRIC'. 'CUBIC-MEAN'
*/
void setSigmaCombiningRule( std::string& sigmaCombiningRule );
void setSigmaCombiningRule( const std::string& sigmaCombiningRule );
/**
* Get sigma combining rule
*
* @return sigmaCombiningRule sigma combining rule: 'ARITHMETIC', 'GEOMETRIC'. 'CUBIC-MEAN'
*/
std::string getSigmaCombiningRule( void ) const;
const std::string& getSigmaCombiningRule( void ) const;
/**
* Set epsilon combining rule
*
* @param epsilonCombiningRule epsilon combining rule: 'ARITHMETIC', 'GEOMETRIC'. 'CUBIC-MEAN'
*/
void setEpsilonCombiningRule( std::string& epsilonCombiningRule );
void setEpsilonCombiningRule( const std::string& epsilonCombiningRule );
/**
* Get epsilon combining rule
*
* @return epsilonCombiningRule epsilon combining rule: 'ARITHMETIC', 'GEOMETRIC'. 'HARMONIC', 'HHG'
*/
std::string getEpsilonCombiningRule( void ) const;
const std::string& getEpsilonCombiningRule( void ) const;
/**
* Set exclusions for specified particle
......
plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp
View file @ 0cae89d3
......@@ -35,6 +35,8 @@
#include "internal/AmoebaMultipoleForceImpl.h"
using namespace OpenMM;
using std::string;
using std::vector;
AmoebaMultipoleForce::AmoebaMultipoleForce() : nonbondedMethod(NoCutoff), pmeBSplineOrder(5), cutoffDistance(1.0), ewaldErrorTol(5e-4), mutualInducedIterationMethod(SOR), mutualInducedMaxIterations(60),
mutualInducedTargetEpsilon(1.0e-05), scalingDistanceCutoff(100.0), electricConstant(138.9354558456) {
......@@ -143,13 +145,13 @@ void AmoebaMultipoleForce::setEwaldErrorTolerance(double tol) {
}
int AmoebaMultipoleForce::addParticle( double charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole, int axisType,
int multipoleAtomId1, int multipoleAtomId2, double thole, double dampingFactor, double polarity) {
multipoles.push_back(MultipoleInfo( charge, molecularDipole, molecularQuadrupole, axisType, multipoleAtomId1, multipoleAtomId2, thole, dampingFactor, polarity));
int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity) {
multipoles.push_back(MultipoleInfo( charge, molecularDipole, molecularQuadrupole, axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY, thole, dampingFactor, polarity));
return multipoles.size()-1;
}
void AmoebaMultipoleForce::getMultipoleParameters(int index, double& charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole,
int& axisType, int& multipoleAtomId1, int& multipoleAtomId2, double& thole, double& dampingFactor, double& polarity ) const {
int& axisType, int& multipoleAtomZ, int& multipoleAtomX, int& multipoleAtomY, double& thole, double& dampingFactor, double& polarity ) const {
charge = multipoles[index].charge;
molecularDipole.resize( 3 );
......@@ -169,8 +171,9 @@ void AmoebaMultipoleForce::getMultipoleParameters(int index, double& charge, std
molecularQuadrupole[8] = multipoles[index].molecularQuadrupole[8];
axisType = multipoles[index].axisType;
multipoleAtomId1 = multipoles[index].multipoleAtomId1;
multipoleAtomId2 = multipoles[index].multipoleAtomId2;
multipoleAtomZ = multipoles[index].multipoleAtomZ;
multipoleAtomX = multipoles[index].multipoleAtomX;
multipoleAtomY = multipoles[index].multipoleAtomY;
thole = multipoles[index].thole;
dampingFactor = multipoles[index].dampingFactor;
......@@ -178,7 +181,7 @@ void AmoebaMultipoleForce::getMultipoleParameters(int index, double& charge, std
}
void AmoebaMultipoleForce::setMultipoleParameters(int index, double charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole,
int axisType, int multipoleAtomId1, int multipoleAtomId2, double thole, double dampingFactor, double polarity ) {
int axisType, int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity ) {
multipoles[index].charge = charge;
......@@ -197,8 +200,9 @@ void AmoebaMultipoleForce::setMultipoleParameters(int index, double charge, std:
multipoles[index].molecularQuadrupole[8] = molecularQuadrupole[8];
multipoles[index].axisType = axisType;
multipoles[index].multipoleAtomId1 = multipoleAtomId1;
multipoles[index].multipoleAtomId2 = multipoleAtomId2;
multipoles[index].multipoleAtomZ = multipoleAtomZ;
multipoles[index].multipoleAtomX = multipoleAtomX;
multipoles[index].multipoleAtomY = multipoleAtomY;
multipoles[index].thole = thole;
multipoles[index].dampingFactor = dampingFactor;
multipoles[index].polarity = polarity;
......
plugins/amoeba/openmmapi/src/AmoebaMultipoleForceImpl.cpp
View file @ 0cae89d3
......@@ -54,12 +54,12 @@ void AmoebaMultipoleForceImpl::initialize(ContextImpl& context) {
for( int ii = 0; ii < system.getNumParticles(); ii++ ){
int index, axisType, multipoleAtomId1, multipoleAtomId2;
int index, axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY;
double charge, thole, dampingFactor, polarity ;
std::vector<double> molecularDipole;
std::vector<double> molecularQuadrupole;
owner.getMultipoleParameters( ii, charge, molecularDipole, molecularQuadrupole, axisType, multipoleAtomId1, multipoleAtomId2,
owner.getMultipoleParameters( ii, charge, molecularDipole, molecularQuadrupole, axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY,
thole, dampingFactor, polarity );
// only 'Z-then-X' or 'Bisector' currently handled
......
plugins/amoeba/openmmapi/src/AmoebaTorsionForce.cpp
View file @ 0cae89d3
......@@ -40,7 +40,7 @@ AmoebaTorsionForce::AmoebaTorsionForce() {
}
int AmoebaTorsionForce::addTorsion(int particle1, int particle2, int particle3, int particle4,
std::vector<double> torsion1, std::vector<double> torsion2, std::vector<double> torsion3 ) {
std::vector<double>& torsion1, std::vector<double>& torsion2, std::vector<double>& torsion3 ) {
torsions.push_back(TorsionInfo(particle1, particle2, particle3, particle4, torsion1, torsion2, torsion3));
return torsions.size()-1;
}
......@@ -64,7 +64,7 @@ void AmoebaTorsionForce::getTorsionParameters(int index, int& particle1, int& pa
}
void AmoebaTorsionForce::setTorsionParameters(int index, int particle1, int particle2, int particle3, int particle4,
std::vector<double> torsion1, std::vector<double> torsion2, std::vector<double> torsion3 ) {
std::vector<double>& torsion1, std::vector<double>& torsion2, std::vector<double>& torsion3 ) {
torsions[index].particle1 = particle1;
torsions[index].particle2 = particle2;
......
plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp
View file @ 0cae89d3
......@@ -35,6 +35,8 @@
#include "internal/AmoebaVdwForceImpl.h"
using namespace OpenMM;
using std::string;
using std::vector;
AmoebaVdwForce::AmoebaVdwForce() : usePBC(0), cutoff(1.0e+10) {
}
......@@ -62,19 +64,19 @@ void AmoebaVdwForce::setParticleParameters(int particleIndex, int ivIndex, int c
parameters[particleIndex].reductionFactor = reductionFactor;
}
void AmoebaVdwForce::setSigmaCombiningRule(std::string& inputSigmaCombiningRule ) {
void AmoebaVdwForce::setSigmaCombiningRule( const std::string& inputSigmaCombiningRule ) {
sigmaCombiningRule = inputSigmaCombiningRule;
}
std::string AmoebaVdwForce::getSigmaCombiningRule( void ) const {
const std::string& AmoebaVdwForce::getSigmaCombiningRule( void ) const {
return sigmaCombiningRule;
}
void AmoebaVdwForce::setEpsilonCombiningRule(std::string& inputEpsilonCombiningRule ) {
void AmoebaVdwForce::setEpsilonCombiningRule( const std::string& inputEpsilonCombiningRule ) {
epsilonCombiningRule = inputEpsilonCombiningRule;
}
std::string AmoebaVdwForce::getEpsilonCombiningRule( void ) const {
const std::string& AmoebaVdwForce::getEpsilonCombiningRule( void ) const {
return epsilonCombiningRule;
}
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernelFactory.cpp
View file @ 0cae89d3
......@@ -38,7 +38,9 @@ extern "C" void OPENMMCUDA_EXPORT registerKernelFactories() {
for( int ii = 0; ii < Platform::getNumPlatforms(); ii++ ){
Platform& platform = Platform::getPlatform(ii);
if( platform.getName() == "Cuda" ){
AmoebaCudaKernelFactory* factory = new AmoebaCudaKernelFactory();
platform.registerKernelFactory(CalcAmoebaHarmonicBondForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaHarmonicAngleForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaHarmonicInPlaneAngleForceKernel::Name(), factory);
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 0cae89d3
......@@ -44,7 +44,9 @@ extern "C" int gpuSetConstants( gpuContext gpu );
using namespace OpenMM;
using namespace std;
// ***************************************************************************
/* -------------------------------------------------------------------------- *
* Calculates bonded forces *
* -------------------------------------------------------------------------- */
static void computeAmoebaLocalForces( AmoebaCudaData& data ) {
amoebaGpuContext gpu = data.getAmoebaGpu();
......@@ -58,6 +60,10 @@ static void computeAmoebaLocalForces( AmoebaCudaData& data ) {
}
/* -------------------------------------------------------------------------- *
* AmoebaHarmonicBond *
* -------------------------------------------------------------------------- */
class CudaCalcAmoebaHarmonicBondForceKernel::ForceInfo : public CudaForceInfo {
public:
ForceInfo(const AmoebaHarmonicBondForce& force) : force(force) {
......@@ -127,6 +133,10 @@ double CudaCalcAmoebaHarmonicBondForceKernel::execute(ContextImpl& context, bool
return 0.0;
}
/* -------------------------------------------------------------------------- *
* AmoebaHarmonicAngle *
* -------------------------------------------------------------------------- */
class CudaCalcAmoebaHarmonicAngleForceKernel::ForceInfo : public CudaForceInfo {
public:
ForceInfo(const AmoebaHarmonicAngleForce& force) : force(force) {
......@@ -195,6 +205,10 @@ double CudaCalcAmoebaHarmonicAngleForceKernel::execute(ContextImpl& context, boo
return 0.0;
}
/* -------------------------------------------------------------------------- *
* AmoebaHarmonicInPlaneAngle *
* -------------------------------------------------------------------------- */
class CudaCalcAmoebaHarmonicInPlaneAngleForceKernel::ForceInfo : public CudaForceInfo {
public:
ForceInfo(const AmoebaHarmonicInPlaneAngleForce& force) : force(force) {
......@@ -267,6 +281,10 @@ double CudaCalcAmoebaHarmonicInPlaneAngleForceKernel::execute(ContextImpl& conte
return 0.0;
}
/* -------------------------------------------------------------------------- *
* AmoebaHarmonicTorsion *
* -------------------------------------------------------------------------- */
class CudaCalcAmoebaTorsionForceKernel::ForceInfo : public CudaForceInfo {
public:
ForceInfo(const AmoebaTorsionForce& force) : force(force) {
......@@ -358,6 +376,10 @@ double CudaCalcAmoebaTorsionForceKernel::execute(ContextImpl& context, bool incl
return 0.0;
}
/* -------------------------------------------------------------------------- *
* AmoebaHarmonicPiTorsion *
* -------------------------------------------------------------------------- */
class CudaCalcAmoebaPiTorsionForceKernel::ForceInfo : public CudaForceInfo {
public:
ForceInfo(const AmoebaPiTorsionForce& force) : force(force) {
......@@ -429,6 +451,10 @@ double CudaCalcAmoebaPiTorsionForceKernel::execute(ContextImpl& context, bool in
return 0.0;
}
/* -------------------------------------------------------------------------- *
* AmoebaStretchBend *
* -------------------------------------------------------------------------- */
class CudaCalcAmoebaStretchBendForceKernel::ForceInfo : public CudaForceInfo {
public:
ForceInfo(const AmoebaStretchBendForce& force) : force(force) {
......@@ -499,6 +525,10 @@ double CudaCalcAmoebaStretchBendForceKernel::execute(ContextImpl& context, bool
return 0.0;
}
/* -------------------------------------------------------------------------- *
* AmoebaOutOfPlaneBend *
* -------------------------------------------------------------------------- */
class CudaCalcAmoebaOutOfPlaneBendForceKernel::ForceInfo : public CudaForceInfo {
public:
ForceInfo(const AmoebaOutOfPlaneBendForce& force) : force(force) {
......@@ -569,6 +599,10 @@ double CudaCalcAmoebaOutOfPlaneBendForceKernel::execute(ContextImpl& context, bo
return 0.0;
}
/* -------------------------------------------------------------------------- *
* AmoebaTorsionTorsion *
* -------------------------------------------------------------------------- */
class CudaCalcAmoebaTorsionTorsionForceKernel::ForceInfo : public CudaForceInfo {
public:
ForceInfo(const AmoebaTorsionTorsionForce& force) : force(force) {
......@@ -716,10 +750,10 @@ public:
}
bool areParticlesIdentical(int particle1, int particle2) {
double charge1, charge2, thole1, thole2, damping1, damping2, polarity1, polarity2;
int axis1, axis2, multipole11, multipole12, multipole21, multipole22;
int axis1, axis2, multipole11, multipole12, multipole21, multipole22, multipole31, multipole32;
vector<double> dipole1, dipole2, quadrupole1, quadrupole2;
force.getMultipoleParameters(particle1, charge1, dipole1, quadrupole1, axis1, multipole11, multipole21, thole1, damping1, polarity1);
force.getMultipoleParameters(particle2, charge2, dipole2, quadrupole2, axis2, multipole12, multipole22, thole2, damping2, polarity2);
force.getMultipoleParameters(particle1, charge1, dipole1, quadrupole1, axis1, multipole11, multipole21, multipole31, thole1, damping1, polarity1);
force.getMultipoleParameters(particle2, charge2, dipole2, quadrupole2, axis2, multipole12, multipole22, multipole32, thole2, damping2, polarity2);
if (charge1 != charge2 || thole1 != thole2 || damping1 != damping2 || polarity1 != polarity2 || axis1 != axis2)
return false;
for (int i = 0; i < (int) dipole1.size(); ++i)
......@@ -756,8 +790,9 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
std::vector<float> dampingFactors(numMultipoles);
std::vector<float> polarity(numMultipoles);
std::vector<int> axisTypes(numMultipoles);
std::vector<int> multipoleAtomId1s(numMultipoles);
std::vector<int> multipoleAtomId2s(numMultipoles);
std::vector<int> multipoleAtomZs(numMultipoles);
std::vector<int> multipoleAtomXs(numMultipoles);
std::vector<int> multipoleAtomYs(numMultipoles);
std::vector< std::vector< std::vector<int> > > multipoleAtomCovalentInfo(numMultipoles);
std::vector<int> minCovalentIndices(numMultipoles);
std::vector<int> minCovalentPolarizationIndices(numMultipoles);
......@@ -786,17 +821,18 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
// multipoles
int axisType, multipoleAtomId1, multipoleAtomId2;
int axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY;
double charge, tholeD, dampingFactorD, polarityD;
std::vector<double> dipolesD;
std::vector<double> quadrupolesD;
force.getMultipoleParameters(i, charge, dipolesD, quadrupolesD, axisType, multipoleAtomId1, multipoleAtomId2,
force.getMultipoleParameters(i, charge, dipolesD, quadrupolesD, axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY,
tholeD, dampingFactorD, polarityD );
totalCharge += charge;
axisTypes[i] = axisType;
multipoleAtomId1s[i] = multipoleAtomId1;
multipoleAtomId2s[i] = multipoleAtomId2;
multipoleAtomZs[i] = multipoleAtomZ;
multipoleAtomXs[i] = multipoleAtomX;
multipoleAtomYs[i] = multipoleAtomY;
charges[i] = static_cast<float>(charge);
tholes[i] = static_cast<float>(tholeD);
......@@ -847,7 +883,7 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
throw OpenMMException("AmoebaMultipoleForce nonbonded method not recognized.\n");
}
gpuSetAmoebaMultipoleParameters(data.getAmoebaGpu(), charges, dipoles, quadrupoles, axisTypes, multipoleAtomId1s, multipoleAtomId2s,
gpuSetAmoebaMultipoleParameters(data.getAmoebaGpu(), charges, dipoles, quadrupoles, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
tholes, scalingDistanceCutoff, dampingFactors, polarity,
multipoleAtomCovalentInfo, covalentDegree, minCovalentIndices, minCovalentPolarizationIndices, (maxCovalentRange+2),
static_cast<int>(force.getMutualInducedIterationMethod()),
......@@ -957,6 +993,10 @@ static void computeAmoebaVdwForce( AmoebaCudaData& data ) {
kCalculateAmoebaVdw14_7Forces(gpu, data.getApplyCutoff());
}
/* -------------------------------------------------------------------------- *
* AmoebaVdw *
* -------------------------------------------------------------------------- */
class CudaCalcAmoebaVdwForceKernel::ForceInfo : public CudaForceInfo {
public:
ForceInfo(const AmoebaVdwForce& force) : force(force) {
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaCudaGpu.cpp
View file @ 0cae89d3
......@@ -75,14 +75,14 @@ amoebaGpuContext amoebaGpuInit( _gpuContext* gpu )
amoebaGpu->amoebaSim.numberOfAtoms = gpu->natoms;
amoebaGpu->amoebaSim.paddedNumberOfAtoms = gpu->sim.paddedNumberOfAtoms;
amoebaGpu->psThetai1 = NULL;
amoebaGpu->psThetai2 = NULL;
amoebaGpu->psThetai3 = NULL;
amoebaGpu->psIgrid = NULL;
amoebaGpu->psPhi = NULL;
amoebaGpu->psPhid = NULL;
amoebaGpu->psPhip = NULL;
amoebaGpu->psPhidp = NULL;
amoebaGpu->psThetai1 = NULL;
amoebaGpu->psThetai2 = NULL;
amoebaGpu->psThetai3 = NULL;
amoebaGpu->psIgrid = NULL;
amoebaGpu->psPhi = NULL;
amoebaGpu->psPhid = NULL;
amoebaGpu->psPhip = NULL;
amoebaGpu->psPhidp = NULL;
return amoebaGpu;
}
......@@ -1446,7 +1446,7 @@ void gpuKirkwoodAllocate( amoebaGpuContext amoebaGpu )
extern "C"
void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vector<float>& charges, const std::vector<float>& dipoles, const std::vector<float>& quadrupoles,
const std::vector<int>& axisType, const std::vector<int>& multipoleParticleId1, const std::vector<int>& multipoleParticleId2,
const std::vector<int>& axisType, const std::vector<int>& multipoleParticleZ, const std::vector<int>& multipoleParticleX, const std::vector<int>& multipoleParticleY,
const std::vector<float>& tholes, float scalingDistanceCutoff,const std::vector<float>& dampingFactors, const std::vector<float>& polarity,
const std::vector< std::vector< std::vector<int> > >& multipoleParticleCovalentInfo, const std::vector<int>& covalentDegree,
const std::vector<int>& minCovalentIndices, const std::vector<int>& minCovalentPolarizationIndices, int maxCovalentRange,
......@@ -1562,10 +1562,10 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
// axis type & multipole particles ids
amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].x = multipoleParticleId1[ii];
amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].y = multipoleParticleId2[ii];
amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].x = multipoleParticleZ[ii];
amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].y = multipoleParticleX[ii];
amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][ii].w = axisType[ii];
int axisParticleIndex = multipoleParticleId1[ii];
int axisParticleIndex = multipoleParticleZ[ii];
if( maxIndices[axisParticleIndex] < ii ){
maxIndices[axisParticleIndex] = ii;
}
......@@ -1573,7 +1573,7 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][axisParticleIndex].z = ii;
}
axisParticleIndex = multipoleParticleId2[ii];
axisParticleIndex = multipoleParticleX[ii];
if( maxIndices[axisParticleIndex] < ii ){
maxIndices[axisParticleIndex] = ii;
}
......@@ -1583,11 +1583,11 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
if( 0 && amoebaGpu->log )
fprintf( amoebaGpu->log, "Z1 %4d [%4d %4d] %4d %4d %4d %4d %d %d\n", ii,
multipoleParticleId1[ii], multipoleParticleId2[ii],
amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][multipoleParticleId1[ii]].z,
maxIndices[multipoleParticleId1[ii]],
amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][multipoleParticleId2[ii]].z,
maxIndices[multipoleParticleId2[ii]],
multipoleParticleZ[ii], multipoleParticleX[ii],
amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][multipoleParticleZ[ii]].z,
maxIndices[multipoleParticleZ[ii]],
amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][multipoleParticleX[ii]].z,
maxIndices[multipoleParticleX[ii]],
amoebaGpu->psMultipoleParticlesIdsAndAxisType->_pSysStream[0][0].z, maxIndices[0] );
// charges
......@@ -1991,7 +1991,6 @@ void gpuSetAmoebaVdwParameters( amoebaGpuContext amoebaGpu,
// ---------------------------------------------------------------------------------------
gpuContext gpu = amoebaGpu->gpuContext;
amoebaGpu->paddedNumberOfAtoms = amoebaGpu->gpuContext->sim.paddedNumberOfAtoms;
unsigned int particles = sigmas.size();
......@@ -2162,7 +2161,11 @@ void gpuSetAmoebaVdwParameters( amoebaGpuContext amoebaGpu,
amoebaGpu->vdwSigmaCombiningRule, amoebaGpu->vdwEpsilonCombiningRule);
for (unsigned int ii = 0; ii < gpu->natoms; ii++)
{
(void) fprintf( amoebaGpu->log, "%5u %15.7e %15.7e\n", ii, sigmas[ii], epsilons[ii] );
(void) fprintf( amoebaGpu->log, "%5u %15.7e %15.7e Ex[", ii, sigmas[ii], epsilons[ii] );
for (unsigned int jj = 0; jj < allExclusions[ii].size(); jj++){
(void) fprintf( amoebaGpu->log, "%d ", allExclusions[ii][jj] );
}
(void) fprintf( amoebaGpu->log, "]\n", ii, sigmas[ii], epsilons[ii] );
if( ii == maxPrint && ii < (amoebaGpu->paddedNumberOfAtoms - maxPrint) )
{
ii = (amoebaGpu->paddedNumberOfAtoms - maxPrint);
......@@ -2172,6 +2175,7 @@ void gpuSetAmoebaVdwParameters( amoebaGpuContext amoebaGpu,
}
#endif
#undef AMOEBA_DEBUG
}
extern "C"
......@@ -2534,7 +2538,7 @@ extern "C"
void gpuSetAmoebaWcaDispersionParameters( amoebaGpuContext amoebaGpu,
const std::vector<float>& radii,
const std::vector<float>& epsilons,
const float totalMaxWcaDisperionEnergy,
const float totalMaxWcaDispersionEnergy,
const float epso, const float epsh, const float rmino, const float rminh,
const float awater, const float shctd, const float dispoff )
{
......@@ -2563,7 +2567,7 @@ void gpuSetAmoebaWcaDispersionParameters( amoebaGpuContext amoebaGpu,
amoebaGpu->psWcaDispersionRadiusEpsilon->_pSysStream[0][ii].y = 0.0f;
}
amoebaGpu->psWcaDispersionRadiusEpsilon->Upload();
amoebaGpu->amoebaSim.totalMaxWcaDispersionEnergy = totalMaxWcaDisperionEnergy;
amoebaGpu->amoebaSim.totalMaxWcaDispersionEnergy = totalMaxWcaDispersionEnergy;
amoebaGpu->amoebaSim.epso = epso;
amoebaGpu->amoebaSim.epsh = epsh;
amoebaGpu->amoebaSim.rmino = rmino;
......@@ -2589,6 +2593,9 @@ void gpuSetAmoebaWcaDispersionParameters( amoebaGpuContext amoebaGpu,
}
#endif
amoebaGpuBuildOutputBuffers( amoebaGpu );
amoebaGpuBuildThreadBlockWorkList( amoebaGpu );
}
extern "C"
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaGpuTypes.h
View file @ 0cae89d3
......@@ -294,7 +294,7 @@ void gpuSetAmoebaTorsionTorsionGrids(amoebaGpuContext gpu, const std::vector< st
extern "C"
void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vector<float>& charges, const std::vector<float>& dipoles, const std::vector<float>& quadrupoles,
const std::vector<int>& axisType, const std::vector<int>& multipoleAtomId1, const std::vector<int>& multipoleAtomId2,
const std::vector<int>& axisType, const std::vector<int>& multipoleAtomZ, const std::vector<int>& multipoleAtomX, const std::vector<int>& multipoleAtomY,
const std::vector<float>& tholes, float scalingDistanceCutoff,const std::vector<float>& dampingFactors, const std::vector<float>& polarity,
const std::vector< std::vector< std::vector<int> > >& multipoleAtomCovalentInfo, const std::vector<int>& covalentDegree,
const std::vector<int>& minCovalentIndices, const std::vector<int>& minCovalentPolarizationIndices, int maxCovalentRange,
......
plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaPmeDirectElectrostatic.cu
View file @ 0cae89d3
......@@ -177,23 +177,21 @@ __device__ void calculatePmeDirectElectrostaticPairIxn_kernel( PmeDirectElectros
// set the permanent multipole and induced dipole values;
float pdi = atomI.damp;
float pti = atomI.thole;
float ci = atomI.q;
float di1 = atomI.labFrameDipole[0];
float di2 = atomI.labFrameDipole[1];
float di3 = atomI.labFrameDipole[2];
float di1 = atomI.labFrameDipole[0];
float di2 = atomI.labFrameDipole[1];
float di3 = atomI.labFrameDipole[2];
float qi1 = atomI.labFrameQuadrupole[0];
float qi2 = atomI.labFrameQuadrupole[1];
float qi3 = atomI.labFrameQuadrupole[2];
float qi4 = atomI.labFrameQuadrupole[3];
float qi5 = atomI.labFrameQuadrupole[4];
float qi6 = atomI.labFrameQuadrupole[5];
float qi7 = atomI.labFrameQuadrupole[6];
float qi8 = atomI.labFrameQuadrupole[7];
float qi9 = atomI.labFrameQuadrupole[8];
float qi1 = atomI.labFrameQuadrupole[0];
float qi2 = atomI.labFrameQuadrupole[1];
float qi3 = atomI.labFrameQuadrupole[2];
float qi4 = atomI.labFrameQuadrupole[3];
float qi5 = atomI.labFrameQuadrupole[4];
float qi6 = atomI.labFrameQuadrupole[5];
float qi7 = atomI.labFrameQuadrupole[6];
float qi8 = atomI.labFrameQuadrupole[7];
float qi9 = atomI.labFrameQuadrupole[8];
float dk1 = atomJ.labFrameDipole[0];
float dk2 = atomJ.labFrameDipole[1];
......@@ -223,19 +221,19 @@ __device__ void calculatePmeDirectElectrostaticPairIxn_kernel( PmeDirectElectros
float exp2a = exp(-(ralpha*ralpha));
alsq2n *= alsq2;
float bn1 = (bn0+alsq2n*exp2a)/r2;
float bn1 = (bn0+alsq2n*exp2a)/r2;
alsq2n *= alsq2;
float bn2 = (3.0f*bn1+alsq2n*exp2a)/r2;
float bn2 = (3.0f*bn1+alsq2n*exp2a)/r2;
alsq2n *= alsq2;
float bn3 = (5.0f*bn2+alsq2n*exp2a)/r2;
float bn3 = (5.0f*bn2+alsq2n*exp2a)/r2;
alsq2n *= alsq2;
float bn4 = (7.0f*bn3+alsq2n*exp2a)/r2;
float bn4 = (7.0f*bn3+alsq2n*exp2a)/r2;
alsq2n *= alsq2;
float bn5 = (9.0f*bn4+alsq2n*exp2a)/r2;
float bn5 = (9.0f*bn4+alsq2n*exp2a)/r2;
// apply Thole polarization damping to scale factors;
......@@ -261,11 +259,9 @@ __device__ void calculatePmeDirectElectrostaticPairIxn_kernel( PmeDirectElectros
float ddsc72 = 0.0f;
float ddsc73 = 0.0f;
float pdk = atomJ.damp;
float ptk = atomJ.thole;
float damp = pdi*pdk;
float damp = atomI.damp*atomJ.damp;
if( damp != 0.0f ){
float pgamma = pti < ptk ? pti : ptk;
float pgamma = atomI.thole < atomJ.thole ? atomI.thole : atomJ.thole;
float ratio = r/damp;
damp = -pgamma * ratio*ratio*ratio;
if( damp > -50.0f ){
......@@ -308,120 +304,136 @@ __device__ void calculatePmeDirectElectrostaticPairIxn_kernel( PmeDirectElectros
float dixdk2 = di3*dk1 - di1*dk3;
float dixdk3 = di1*dk2 - di2*dk1;
float dixuk1 = di2*atomJ.inducedDipole[2] - di3*atomJ.inducedDipole[1];
float dixuk2 = di3*atomJ.inducedDipole[0] - di1*atomJ.inducedDipole[2];
float dixuk3 = di1*atomJ.inducedDipole[1] - di2*atomJ.inducedDipole[0];
float dkxui1 = dk2*atomI.inducedDipole[2] - dk3*atomI.inducedDipole[1];
float dkxui2 = dk3*atomI.inducedDipole[0] - dk1*atomI.inducedDipole[2];
float dkxui3 = dk1*atomI.inducedDipole[1] - dk2*atomI.inducedDipole[0];
float dixuk1 = di2*atomJ.inducedDipole[2] - di3*atomJ.inducedDipole[1];
float dixuk2 = di3*atomJ.inducedDipole[0] - di1*atomJ.inducedDipole[2];
float dixuk3 = di1*atomJ.inducedDipole[1] - di2*atomJ.inducedDipole[0];
float dkxui1 = dk2*atomI.inducedDipole[2] - dk3*atomI.inducedDipole[1];
float dkxui2 = dk3*atomI.inducedDipole[0] - dk1*atomI.inducedDipole[2];
float dkxui3 = dk1*atomI.inducedDipole[1] - dk2*atomI.inducedDipole[0];
float dixukp1 = di2*atomJ.inducedDipoleP[2] - di3*atomJ.inducedDipoleP[1];
float dixukp2 = di3*atomJ.inducedDipoleP[0] - di1*atomJ.inducedDipoleP[2];
float dixukp3 = di1*atomJ.inducedDipoleP[1] - di2*atomJ.inducedDipoleP[0];
float dkxuip1 = dk2*atomI.inducedDipoleP[2] - dk3*atomI.inducedDipoleP[1];
float dkxuip2 = dk3*atomI.inducedDipoleP[0] - dk1*atomI.inducedDipoleP[2];
float dkxuip3 = dk1*atomI.inducedDipoleP[1] - dk2*atomI.inducedDipoleP[0];
float dixr1 = di2*zr - di3*yr;
float dixr2 = di3*xr - di1*zr;
float dixr3 = di1*yr - di2*xr;
float dkxr1 = dk2*zr - dk3*yr;
float dkxr2 = dk3*xr - dk1*zr;
float dkxr3 = dk1*yr - dk2*xr;
float qir1 = qi1*xr + qi4*yr + qi7*zr;
float qir2 = qi2*xr + qi5*yr + qi8*zr;
float qir3 = qi3*xr + qi6*yr + qi9*zr;
float qkr1 = qk1*xr + qk4*yr + qk7*zr;
float qkr2 = qk2*xr + qk5*yr + qk8*zr;
float qkr3 = qk3*xr + qk6*yr + qk9*zr;
float qiqkr1 = qi1*qkr1 + qi4*qkr2 + qi7*qkr3;
float qiqkr2 = qi2*qkr1 + qi5*qkr2 + qi8*qkr3;
float qiqkr3 = qi3*qkr1 + qi6*qkr2 + qi9*qkr3;
float qkqir1 = qk1*qir1 + qk4*qir2 + qk7*qir3;
float qkqir2 = qk2*qir1 + qk5*qir2 + qk8*qir3;
float qkqir3 = qk3*qir1 + qk6*qir2 + qk9*qir3;
float qixqk1 = qi2*qk3 + qi5*qk6 + qi8*qk9
- qi3*qk2 - qi6*qk5 - qi9*qk8;
float qixqk2 = qi3*qk1 + qi6*qk4 + qi9*qk7
- qi1*qk3 - qi4*qk6 - qi7*qk9;
float qixqk3 = qi1*qk2 + qi4*qk5 + qi7*qk8
- qi2*qk1 - qi5*qk4 - qi8*qk7;
float qixqk1 = qi2*qk3 + qi5*qk6 + qi8*qk9 - qi3*qk2 - qi6*qk5 - qi9*qk8;
float qixqk2 = qi3*qk1 + qi6*qk4 + qi9*qk7 - qi1*qk3 - qi4*qk6 - qi7*qk9;
float qixqk3 = qi1*qk2 + qi4*qk5 + qi7*qk8 - qi2*qk1 - qi5*qk4 - qi8*qk7;
float rxqir1 = yr*qir3 - zr*qir2;
float rxqir2 = zr*qir1 - xr*qir3;
float rxqir3 = xr*qir2 - yr*qir1;
float rxqkr1 = yr*qkr3 - zr*qkr2;
float rxqkr2 = zr*qkr1 - xr*qkr3;
float rxqkr3 = xr*qkr2 - yr*qkr1;
float rxqikr1 = yr*qiqkr3 - zr*qiqkr2;
float rxqikr2 = zr*qiqkr1 - xr*qiqkr3;
float rxqikr3 = xr*qiqkr2 - yr*qiqkr1;
float rxqkir1 = yr*qkqir3 - zr*qkqir2;
float rxqkir2 = zr*qkqir1 - xr*qkqir3;
float rxqkir3 = xr*qkqir2 - yr*qkqir1;
float qkrxqir1 = qkr2*qir3 - qkr3*qir2;
float qkrxqir2 = qkr3*qir1 - qkr1*qir3;
float qkrxqir3 = qkr1*qir2 - qkr2*qir1;
float qidk1 = qi1*dk1 + qi4*dk2 + qi7*dk3;
float qidk2 = qi2*dk1 + qi5*dk2 + qi8*dk3;
float qidk3 = qi3*dk1 + qi6*dk2 + qi9*dk3;
float qkdi1 = qk1*di1 + qk4*di2 + qk7*di3;
float qkdi2 = qk2*di1 + qk5*di2 + qk8*di3;
float qkdi3 = qk3*di1 + qk6*di2 + qk9*di3;
float qiuk1 = qi1*atomJ.inducedDipole[0] + qi4*atomJ.inducedDipole[1]
+ qi7*atomJ.inducedDipole[2];
float qiuk2 = qi2*atomJ.inducedDipole[0] + qi5*atomJ.inducedDipole[1]
+ qi8*atomJ.inducedDipole[2];
float qiuk3 = qi3*atomJ.inducedDipole[0] + qi6*atomJ.inducedDipole[1]
+ qi9*atomJ.inducedDipole[2];
float qkui1 = qk1*atomI.inducedDipole[0] + qk4*atomI.inducedDipole[1]
+ qk7*atomI.inducedDipole[2];
float qkui2 = qk2*atomI.inducedDipole[0] + qk5*atomI.inducedDipole[1]
+ qk8*atomI.inducedDipole[2];
float qkui3 = qk3*atomI.inducedDipole[0] + qk6*atomI.inducedDipole[1]
+ qk9*atomI.inducedDipole[2];
float qiukp1 = qi1*atomJ.inducedDipoleP[0] + qi4*atomJ.inducedDipoleP[1]
+ qi7*atomJ.inducedDipoleP[2];
float qiukp2 = qi2*atomJ.inducedDipoleP[0] + qi5*atomJ.inducedDipoleP[1]
+ qi8*atomJ.inducedDipoleP[2];
float qiukp3 = qi3*atomJ.inducedDipoleP[0] + qi6*atomJ.inducedDipoleP[1]
+ qi9*atomJ.inducedDipoleP[2];
float qkuip1 = qk1*atomI.inducedDipoleP[0] + qk4*atomI.inducedDipoleP[1]
+ qk7*atomI.inducedDipoleP[2];
float qkuip2 = qk2*atomI.inducedDipoleP[0] + qk5*atomI.inducedDipoleP[1]
+ qk8*atomI.inducedDipoleP[2];
float qkuip3 = qk3*atomI.inducedDipoleP[0] + qk6*atomI.inducedDipoleP[1]
+ qk9*atomI.inducedDipoleP[2];
float qiuk1 = qi1*atomJ.inducedDipole[0] + qi4*atomJ.inducedDipole[1] + qi7*atomJ.inducedDipole[2];
float qiuk2 = qi2*atomJ.inducedDipole[0] + qi5*atomJ.inducedDipole[1] + qi8*atomJ.inducedDipole[2];
float qiuk3 = qi3*atomJ.inducedDipole[0] + qi6*atomJ.inducedDipole[1] + qi9*atomJ.inducedDipole[2];
float qkui1 = qk1*atomI.inducedDipole[0] + qk4*atomI.inducedDipole[1] + qk7*atomI.inducedDipole[2];
float qkui2 = qk2*atomI.inducedDipole[0] + qk5*atomI.inducedDipole[1] + qk8*atomI.inducedDipole[2];
float qkui3 = qk3*atomI.inducedDipole[0] + qk6*atomI.inducedDipole[1] + qk9*atomI.inducedDipole[2];
float qiukp1 = qi1*atomJ.inducedDipoleP[0] + qi4*atomJ.inducedDipoleP[1] + qi7*atomJ.inducedDipoleP[2];
float qiukp2 = qi2*atomJ.inducedDipoleP[0] + qi5*atomJ.inducedDipoleP[1] + qi8*atomJ.inducedDipoleP[2];
float qiukp3 = qi3*atomJ.inducedDipoleP[0] + qi6*atomJ.inducedDipoleP[1] + qi9*atomJ.inducedDipoleP[2];
float qkuip1 = qk1*atomI.inducedDipoleP[0] + qk4*atomI.inducedDipoleP[1] + qk7*atomI.inducedDipoleP[2];
float qkuip2 = qk2*atomI.inducedDipoleP[0] + qk5*atomI.inducedDipoleP[1] + qk8*atomI.inducedDipoleP[2];
float qkuip3 = qk3*atomI.inducedDipoleP[0] + qk6*atomI.inducedDipoleP[1] + qk9*atomI.inducedDipoleP[2];
float dixqkr1 = di2*qkr3 - di3*qkr2;
float dixqkr2 = di3*qkr1 - di1*qkr3;
float dixqkr3 = di1*qkr2 - di2*qkr1;
float dkxqir1 = dk2*qir3 - dk3*qir2;
float dkxqir2 = dk3*qir1 - dk1*qir3;
float dkxqir3 = dk1*qir2 - dk2*qir1;
float uixqkr1 = atomI.inducedDipole[1]*qkr3 - atomI.inducedDipole[2]*qkr2;
float uixqkr2 = atomI.inducedDipole[2]*qkr1 - atomI.inducedDipole[0]*qkr3;
float uixqkr3 = atomI.inducedDipole[0]*qkr2 - atomI.inducedDipole[1]*qkr1;
float ukxqir1 = atomJ.inducedDipole[1]*qir3 - atomJ.inducedDipole[2]*qir2;
float ukxqir2 = atomJ.inducedDipole[2]*qir1 - atomJ.inducedDipole[0]*qir3;
float ukxqir3 = atomJ.inducedDipole[0]*qir2 - atomJ.inducedDipole[1]*qir1;
float uixqkrp1 = atomI.inducedDipoleP[1]*qkr3 - atomI.inducedDipoleP[2]*qkr2;
float uixqkrp2 = atomI.inducedDipoleP[2]*qkr1 - atomI.inducedDipoleP[0]*qkr3;
float uixqkrp3 = atomI.inducedDipoleP[0]*qkr2 - atomI.inducedDipoleP[1]*qkr1;
float ukxqirp1 = atomJ.inducedDipoleP[1]*qir3 - atomJ.inducedDipoleP[2]*qir2;
float ukxqirp2 = atomJ.inducedDipoleP[2]*qir1 - atomJ.inducedDipoleP[0]*qir3;
float ukxqirp3 = atomJ.inducedDipoleP[0]*qir2 - atomJ.inducedDipoleP[1]*qir1;
float rxqidk1 = yr*qidk3 - zr*qidk2;
float rxqidk2 = zr*qidk1 - xr*qidk3;
float rxqidk3 = xr*qidk2 - yr*qidk1;
float rxqkdi1 = yr*qkdi3 - zr*qkdi2;
float rxqkdi2 = zr*qkdi1 - xr*qkdi3;
float rxqkdi3 = xr*qkdi2 - yr*qkdi1;
float rxqiuk1 = yr*qiuk3 - zr*qiuk2;
float rxqiuk2 = zr*qiuk1 - xr*qiuk3;
float rxqiuk3 = xr*qiuk2 - yr*qiuk1;
float rxqkui1 = yr*qkui3 - zr*qkui2;
float rxqkui2 = zr*qkui1 - xr*qkui3;
float rxqkui3 = xr*qkui2 - yr*qkui1;
float rxqiukp1 = yr*qiukp3 - zr*qiukp2;
float rxqiukp2 = zr*qiukp1 - xr*qiukp3;
float rxqiukp3 = xr*qiukp2 - yr*qiukp1;
float rxqkuip1 = yr*qkuip3 - zr*qkuip2;
float rxqkuip2 = zr*qkuip1 - xr*qkuip3;
float rxqkuip3 = xr*qkuip2 - yr*qkuip1;
......@@ -437,33 +449,29 @@ __device__ void calculatePmeDirectElectrostaticPairIxn_kernel( PmeDirectElectros
float sc8 = qkr1*di1 + qkr2*di2 + qkr3*di3;
float sc9 = qir1*qkr1 + qir2*qkr2 + qir3*qkr3;
float sc10 = qi1*qk1 + qi2*qk2 + qi3*qk3
+ qi4*qk4 + qi5*qk5 + qi6*qk6
+ qi7*qk7 + qi8*qk8 + qi9*qk9;
+ qi4*qk4 + qi5*qk5 + qi6*qk6
+ qi7*qk7 + qi8*qk8 + qi9*qk9;
// calculate the scalar products for induced components
float sci1 = atomI.inducedDipole[0]*dk1 + atomI.inducedDipole[1]*dk2
+ atomI.inducedDipole[2]*dk3 + di1*atomJ.inducedDipole[0]
+ di2*atomJ.inducedDipole[1] + di3*atomJ.inducedDipole[2];
float sci3 = atomI.inducedDipole[0]*xr + atomI.inducedDipole[1]*yr + atomI.inducedDipole[2]*zr;
float sci4 = atomJ.inducedDipole[0]*xr + atomJ.inducedDipole[1]*yr + atomJ.inducedDipole[2]*zr;
float sci7 = qir1*atomJ.inducedDipole[0] + qir2*atomJ.inducedDipole[1]
+ qir3*atomJ.inducedDipole[2];
float sci8 = qkr1*atomI.inducedDipole[0] + qkr2*atomI.inducedDipole[1]
+ qkr3*atomI.inducedDipole[2];
float scip1 = atomI.inducedDipoleP[0]*dk1 + atomI.inducedDipoleP[1]*dk2
+ atomI.inducedDipoleP[2]*dk3 + di1*atomJ.inducedDipoleP[0]
+ di2*atomJ.inducedDipoleP[1] + di3*atomJ.inducedDipoleP[2];
float scip2 = atomI.inducedDipole[0]*atomJ.inducedDipoleP[0]+atomI.inducedDipole[1]*atomJ.inducedDipoleP[1]
+ atomI.inducedDipole[2]*atomJ.inducedDipoleP[2]+atomI.inducedDipoleP[0]*atomJ.inducedDipole[0]
+ atomI.inducedDipoleP[1]*atomJ.inducedDipole[1]+atomI.inducedDipoleP[2]*atomJ.inducedDipole[2];
float scip3 = atomI.inducedDipoleP[0]*xr + atomI.inducedDipoleP[1]*yr + atomI.inducedDipoleP[2]*zr;
float scip4 = atomJ.inducedDipoleP[0]*xr + atomJ.inducedDipoleP[1]*yr + atomJ.inducedDipoleP[2]*zr;
float scip7 = qir1*atomJ.inducedDipoleP[0] + qir2*atomJ.inducedDipoleP[1]
+ qir3*atomJ.inducedDipoleP[2];
float scip8 = qkr1*atomI.inducedDipoleP[0] + qkr2*atomI.inducedDipoleP[1]
+ qkr3*atomI.inducedDipoleP[2];
float sci1 = atomI.inducedDipole[0]*dk1 + atomI.inducedDipole[1]*dk2
+ atomI.inducedDipole[2]*dk3 + di1*atomJ.inducedDipole[0]
+ di2*atomJ.inducedDipole[1] + di3*atomJ.inducedDipole[2];
float sci3 = atomI.inducedDipole[0]*xr + atomI.inducedDipole[1]*yr + atomI.inducedDipole[2]*zr;
float sci4 = atomJ.inducedDipole[0]*xr + atomJ.inducedDipole[1]*yr + atomJ.inducedDipole[2]*zr;
float sci7 = qir1*atomJ.inducedDipole[0] + qir2*atomJ.inducedDipole[1] + qir3*atomJ.inducedDipole[2];
float sci8 = qkr1*atomI.inducedDipole[0] + qkr2*atomI.inducedDipole[1] + qkr3*atomI.inducedDipole[2];
float scip1 = atomI.inducedDipoleP[0]*dk1 + atomI.inducedDipoleP[1]*dk2 + atomI.inducedDipoleP[2]*dk3 + di1*atomJ.inducedDipoleP[0] + di2*atomJ.inducedDipoleP[1] + di3*atomJ.inducedDipoleP[2];
float scip2 = atomI.inducedDipole[0]*atomJ.inducedDipoleP[0]+atomI.inducedDipole[1]*atomJ.inducedDipoleP[1]
+ atomI.inducedDipole[2]*atomJ.inducedDipoleP[2]+atomI.inducedDipoleP[0]*atomJ.inducedDipole[0]
+ atomI.inducedDipoleP[1]*atomJ.inducedDipole[1]+atomI.inducedDipoleP[2]*atomJ.inducedDipole[2];
float scip3 = atomI.inducedDipoleP[0]*xr + atomI.inducedDipoleP[1]*yr + atomI.inducedDipoleP[2]*zr;
float scip4 = atomJ.inducedDipoleP[0]*xr + atomJ.inducedDipoleP[1]*yr + atomJ.inducedDipoleP[2]*zr;
float scip7 = qir1*atomJ.inducedDipoleP[0] + qir2*atomJ.inducedDipoleP[1] + qir3*atomJ.inducedDipoleP[2];
float scip8 = qkr1*atomI.inducedDipoleP[0] + qkr2*atomI.inducedDipoleP[1] + qkr3*atomI.inducedDipoleP[2];
// calculate the gl functions for permanent components
......@@ -492,24 +500,17 @@ __device__ void calculatePmeDirectElectrostaticPairIxn_kernel( PmeDirectElectros
// compute the energy contributions for this interaction
float e = bn0*gl0 + bn1*(gl1+gl6)
+ bn2*(gl2+gl7+gl8)
+ bn3*(gl3+gl5) + bn4*gl4;
float ei = 0.5f * (bn1*(gli1+gli6)
+ bn2*(gli2+gli7) + bn3*gli3);
float e = bn0*gl0 + bn1*(gl1+gl6) + bn2*(gl2+gl7+gl8) + bn3*(gl3+gl5) + bn4*gl4;
float ei = 0.5f * (bn1*(gli1+gli6) + bn2*(gli2+gli7) + bn3*gli3);
// get the real energy without any screening function
float erl = rr1*gl0 + rr3*(gl1+gl6)
+ rr5*(gl2+gl7+gl8)
+ rr7*(gl3+gl5) + rr9*gl4;
float erli = 0.5f*(rr3*(gli1+gli6)*psc3
+ rr5*(gli2+gli7)*psc5
+ rr7*gli3*psc7);
e = e - (1.0f-scalingFactors[MScaleIndex])*erl;
ei = ei - erli;
float erl = rr1*gl0 + rr3*(gl1+gl6) + rr5*(gl2+gl7+gl8) + rr7*(gl3+gl5) + rr9*gl4;
float erli = 0.5f*(rr3*(gli1+gli6)*psc3 + rr5*(gli2+gli7)*psc5 + rr7*gli3*psc7);
e = e - (1.0f-scalingFactors[MScaleIndex])*erl;
ei = ei - erli;
*energy = -conversionFactor*(e + ei);
*energy = -conversionFactor*(e + ei);
// increment the total intramolecular energy; assumes;
// intramolecular distances are less than half of cell;
......@@ -638,6 +639,7 @@ __device__ void calculatePmeDirectElectrostaticPairIxn_kernel( PmeDirectElectros
+ 0.5f*gfri4*((qkui1-qiuk1)*psc5
+ (qkuip1-qiukp1)*dsc5)
+ gfri5*qir1 + gfri6*qkr1;
float ftm2ri2 = gfri1*yr + 0.5f*
(- rr3*ck*(atomI.inducedDipole[1]*psc3+atomI.inducedDipoleP[1]*dsc3)
+ rr5*sc4*(atomI.inducedDipole[1]*psc5+atomI.inducedDipoleP[1]*dsc5)
......@@ -652,6 +654,7 @@ __device__ void calculatePmeDirectElectrostaticPairIxn_kernel( PmeDirectElectros
+ 0.5f*gfri4*((qkui2-qiuk2)*psc5
+ (qkuip2-qiukp2)*dsc5)
+ gfri5*qir2 + gfri6*qkr2;
float ftm2ri3 = gfri1*zr + 0.5f*
(- rr3*ck*(atomI.inducedDipole[2]*psc3+atomI.inducedDipoleP[2]*dsc3)
+ rr5*sc4*(atomI.inducedDipole[2]*psc5+atomI.inducedDipoleP[2]*dsc5)
......@@ -950,9 +953,9 @@ __device__ void calculatePmeDirectElectrostaticPairIxn_kernel( PmeDirectElectros
} else {
outputForce[0] = 0.0f;
outputForce[1] = 0.0f;
outputForce[2] = 0.0f;
outputForce[0] = 0.0f;
outputForce[1] = 0.0f;
outputForce[2] = 0.0f;
outputTorque[0].x = 0.0f;
outputTorque[0].y = 0.0f;
......@@ -962,7 +965,7 @@ __device__ void calculatePmeDirectElectrostaticPairIxn_kernel( PmeDirectElectros
outputTorque[1].y = 0.0f;
outputTorque[1].z = 0.0f;
*energy = 0.0f;
*energy = 0.0f;
#ifdef AMOEBA_DEBUG
for( int ii = 0; ii < 5; ii++ ){
......
plugins/amoeba/platforms/cuda/tests/AmoebaTinkerParameterFile.cpp
View file @ 0cae89d3
......@@ -1600,7 +1600,6 @@ static int readAmoebaOutOfPlaneBendParameters( FILE* filePtr, MapStringInt& forc
int particle1, particle2, particle3, particle4;
double k;
outOfPlaneBendForce->getOutOfPlaneBendParameters( ii, particle1, particle2, particle3, particle4, k );
//k *= CalToJoule/(AngstromToNm*AngstromToNm);
k *= CalToJoule;
outOfPlaneBendForce->setOutOfPlaneBendParameters( ii, particle1, particle2, particle3, particle4, k );
}
......@@ -2053,6 +2052,12 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
int axisType = atoi( lineTokens[tokenIndex++].c_str() );
int zAxis = atoi( lineTokens[tokenIndex++].c_str() );
int xAxis = atoi( lineTokens[tokenIndex++].c_str() );
int yAxis;
if( version > 2 ){
yAxis = atoi( lineTokens[tokenIndex++].c_str() );
} else {
yAxis = -1;
}
double pdamp = atof( lineTokens[tokenIndex++].c_str() );
double tholeDamp = atof( lineTokens[tokenIndex++].c_str() );
double polarity = atof( lineTokens[tokenIndex++].c_str() );
......@@ -2069,7 +2074,7 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
quadrupole[6] = atof( lineTokens[tokenIndex++].c_str() );
quadrupole[7] = atof( lineTokens[tokenIndex++].c_str() );
quadrupole[8] = atof( lineTokens[tokenIndex++].c_str() );
multipoleForce->addParticle( charge, dipole, quadrupole, axisType, zAxis, xAxis, tholeDamp, pdamp, polarity );
multipoleForce->addParticle( charge, dipole, quadrupole, axisType, zAxis, xAxis, yAxis, tholeDamp, pdamp, polarity );
} else {
(void) fprintf( log, "%s AmoebaMultipoleForce tokens incomplete at line=%d\n", methodName.c_str(), *lineCount );
(void) fflush( log );
......@@ -2160,11 +2165,11 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
for( int ii = 0; ii < multipoleForce->getNumMultipoles(); ii++ ){
int axisType, zAxis, xAxis;
int axisType, zAxis, xAxis, yAxis;
std::vector<double> dipole;
std::vector<double> quadrupole;
double charge, thole, dampingFactor, polarity;
multipoleForce->getMultipoleParameters( ii, charge, dipole, quadrupole, axisType, zAxis, xAxis, thole, dampingFactor, polarity );
multipoleForce->getMultipoleParameters( ii, charge, dipole, quadrupole, axisType, zAxis, xAxis, yAxis, thole, dampingFactor, polarity );
for( unsigned int jj = 0; jj < dipole.size(); jj++ ){
dipole[jj] *= dipoleConversion;
......@@ -2175,7 +2180,7 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
polarity *= polarityConversion;
dampingFactor *= dampingFactorConversion;
multipoleForce->setMultipoleParameters( ii, charge, dipole, quadrupole, axisType, zAxis, xAxis, thole, dampingFactor, polarity );
multipoleForce->setMultipoleParameters( ii, charge, dipole, quadrupole, axisType, zAxis, xAxis, yAxis, thole, dampingFactor, polarity );
}
} else {
......@@ -2227,11 +2232,11 @@ static int readAmoebaMultipoleParameters( FILE* filePtr, int version, MapStringI
(void) fflush( log );
for( unsigned int ii = 0; ii < arraySize; ii++ ){
int axisType, zAxis, xAxis;
int axisType, zAxis, xAxis, yAxis;
std::vector<double> dipole;
std::vector<double> quadrupole;
double charge, thole, dampingFactor, polarity;
multipoleForce->getMultipoleParameters( ii, charge, dipole, quadrupole, axisType, zAxis, xAxis, thole, dampingFactor, polarity );
multipoleForce->getMultipoleParameters( ii, charge, dipole, quadrupole, axisType, zAxis, xAxis, yAxis, thole, dampingFactor, polarity );
(void) fprintf( log, "%8d %8d %8d %8d q %10.4f thl %10.4f pgm %10.4f pol %10.4f d[%10.4f %10.4f %10.4f]\n",
ii, axisType, zAxis, xAxis, charge, thole, dampingFactor, polarity, dipole[0], dipole[1], dipole[2] );
(void) fprintf( log, " q[%10.4f %10.4f %10.4f] [%10.4f %10.4f %10.4f] [%10.4f %10.4f %10.4f]\n",
......@@ -2361,7 +2366,7 @@ static int readAmoebaGeneralizedKirkwoodParameters( FILE* filePtr, MapStringInt&
}
}
// get supplementart fields
// get supplementary fields
isNotEof = 1;
int totalFields = 2;
......@@ -2850,7 +2855,7 @@ static int readAmoebaWcaDispersionParameters( FILE* filePtr, MapStringInt& force
//static const unsigned int maxPrint = MAX_PRINT;
static const unsigned int maxPrint = 15;
unsigned int arraySize = static_cast<unsigned int>(wcaDispersionForce->getNumParticles());
(void) fprintf( log, "%s: %u sample of AmoebaVdwForce parameters in %s units.\n",
(void) fprintf( log, "%s: %u sample of AmoebaWcaForce parameters in %s units.\n",
methodName.c_str(), arraySize, (useOpenMMUnits ? "OpenMM" : "Amoeba") );
(void) fprintf( log, "Eps[%14.7f %14.7f] Rmin[%14.7f %14.7f]\nAwater %14.7f Shctd %14.7f Dispoff %14.7f Slevy %14.7f\n",
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaPME.cpp
View file @ 0cae89d3
......@@ -65,7 +65,7 @@ void testPMEWater() {
quadrupole[8] = 3.622635596e-3;
double damp = 9.707801995e-01*sqrt(0.1);
double polarity = 0.837*0.001;
mp->addParticle(-0.51966, dipole, quadrupole, 1, 1, 2, 0.39, damp, polarity);
mp->addParticle(-0.51966, dipole, quadrupole, 1, 1, 2, -1, 0.39, damp, polarity);
dipole[0] = -2.042094848e-2;
dipole[2] = -3.078753000e-2;
quadrupole[0] = -3.428482490e-3;
......@@ -75,8 +75,8 @@ void testPMEWater() {
quadrupole[8] = 1.345257001e-2;
damp = 8.897068742e-01*sqrt(0.1);
polarity = 0.496*0.001;
mp->addParticle(0.25983, dipole, quadrupole, 0, 0, 2, 0.39, damp, polarity);
mp->addParticle(0.25983, dipole, quadrupole, 0, 0, 1, 0.39, damp, polarity);
mp->addParticle(0.25983, dipole, quadrupole, 0, 0, 2, -1, 0.39, damp, polarity);
mp->addParticle(0.25983, dipole, quadrupole, 0, 0, 1, -1, 0.39, damp, polarity);
mp->setCutoffDistance(1.0);
std::vector<int> intVector;
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ 0cae89d3
......@@ -423,8 +423,9 @@ void ReferenceCalcAmoebaMultipoleForceKernel::initialize(const System& system, c
dampingFactors.resize(numMultipoles);
polarity.resize(numMultipoles);
axisTypes.resize(numMultipoles);
multipoleAtomId1s.resize(numMultipoles);
multipoleAtomId2s.resize(numMultipoles);
multipoleAtomZs.resize(numMultipoles);
multipoleAtomXs.resize(numMultipoles);
multipoleAtomYs.resize(numMultipoles);
multipoleAtomCovalentInfo.resize(numMultipoles);
int dipoleIndex = 0;
......@@ -435,17 +436,18 @@ void ReferenceCalcAmoebaMultipoleForceKernel::initialize(const System& system, c
// multipoles
int axisType, multipoleAtomId1, multipoleAtomId2;
int axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY;
double charge, tholeD, dampingFactorD, polarityD;
std::vector<double> dipolesD;
std::vector<double> quadrupolesD;
force.getMultipoleParameters(ii, charge, dipolesD, quadrupolesD, axisType, multipoleAtomId1, multipoleAtomId2,
force.getMultipoleParameters(ii, charge, dipolesD, quadrupolesD, axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY,
tholeD, dampingFactorD, polarityD );
totalCharge += charge;
axisTypes[ii] = axisType;
multipoleAtomId1s[ii] = multipoleAtomId1;
multipoleAtomId2s[ii] = multipoleAtomId2;
multipoleAtomZs[ii] = multipoleAtomZ;
multipoleAtomXs[ii] = multipoleAtomX;
multipoleAtomYs[ii] = multipoleAtomY;
charges[ii] = static_cast<RealOpenMM>(charge);
tholes[ii] = static_cast<RealOpenMM>(tholeD);
......@@ -494,7 +496,7 @@ double ReferenceCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bo
RealOpenMM energy = amoebaReferenceMultipoleForce.calculateForceAndEnergy( numMultipoles, posData,
charges, dipoles, quadrupoles, tholes,
dampingFactors, polarity, axisTypes,
multipoleAtomId1s, multipoleAtomId2s,
multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
multipoleAtomCovalentInfo, forceData);
return static_cast<double>(energy);
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ 0cae89d3
......@@ -357,8 +357,9 @@ private:
std::vector<RealOpenMM> dampingFactors;
std::vector<RealOpenMM> polarity;
std::vector<int> axisTypes;
std::vector<int> multipoleAtomId1s;
std::vector<int> multipoleAtomId2s;
std::vector<int> multipoleAtomZs;
std::vector<int> multipoleAtomXs;
std::vector<int> multipoleAtomYs;
std::vector< std::vector< std::vector<int> > > multipoleAtomCovalentInfo;
//int iterativeMethod;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp
View file @ 0cae89d3
......@@ -1443,8 +1443,9 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForce( const MultipoleParticleDat
RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostatic( std::vector<MultipoleParticleData>& particleData,
const std::vector<int>& axisTypes,
const std::vector<int>& multipoleAtomId1s,
const std::vector<int>& multipoleAtomId2s,
const std::vector<int>& multipoleAtomZs,
const std::vector<int>& multipoleAtomXs,
const std::vector<int>& multipoleAtomYs,
RealOpenMM** forces ) const {
// ---------------------------------------------------------------------------------------
......@@ -1523,7 +1524,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffElectrostatic( std::v
// map torques to forces
for( unsigned int ii = 0; ii < particleData.size(); ii++ ){
mapTorqueToForce( particleData[ii], particleData[multipoleAtomId1s[ii]], particleData[multipoleAtomId2s[ii]], axisTypes[ii], torques[ii], forces );
mapTorqueToForce( particleData[ii], particleData[multipoleAtomZs[ii]], particleData[multipoleAtomXs[ii]], axisTypes[ii], torques[ii], forces );
}
// diagnostics
......@@ -1555,8 +1556,9 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffForceAndEnergy( unsig
const std::vector<RealOpenMM>& dampingFactors,
const std::vector<RealOpenMM>& polarity,
const std::vector<int>& axisTypes,
const std::vector<int>& multipoleAtomId1s,
const std::vector<int>& multipoleAtomId2s,
const std::vector<int>& multipoleAtomZs,
const std::vector<int>& multipoleAtomXs,
const std::vector<int>& multipoleAtomYs,
RealOpenMM** forces ){
......@@ -1579,8 +1581,8 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffForceAndEnergy( unsig
tholes, dampingFactors, polarity, particleData );
for( unsigned int ii = 0; ii < numParticles; ii++ ){
if( multipoleAtomId1s[ii] >= 0 && multipoleAtomId2s[ii] >= 0 ){
applyRotationMatrix( particleData[ii], particleData[multipoleAtomId1s[ii]], particleData[multipoleAtomId2s[ii]], axisTypes[ii] );
if( multipoleAtomZs[ii] >= 0 && multipoleAtomXs[ii] >= 0 ){
applyRotationMatrix( particleData[ii], particleData[multipoleAtomZs[ii]], particleData[multipoleAtomXs[ii]], axisTypes[ii] );
}
}
......@@ -1637,7 +1639,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateNoCutoffForceAndEnergy( unsig
throw OpenMMException(message.str());
}
RealOpenMM energy = calculateNoCutoffElectrostatic( particleData, axisTypes, multipoleAtomId1s, multipoleAtomId2s, forces );
RealOpenMM energy = calculateNoCutoffElectrostatic( particleData, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, forces );
return energy;
......@@ -1652,8 +1654,9 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateForceAndEnergy( int numPartic
const std::vector<RealOpenMM>& dampingFactors,
const std::vector<RealOpenMM>& polarity,
const std::vector<int>& axisTypes,
const std::vector<int>& multipoleAtomId1s,
const std::vector<int>& multipoleAtomId2s,
const std::vector<int>& multipoleAtomZs,
const std::vector<int>& multipoleAtomXs,
const std::vector<int>& multipoleAtomYs,
const std::vector< std::vector< std::vector<int> > >& multipoleParticleCovalentInfo,
RealOpenMM** forces ){
......@@ -1662,7 +1665,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateForceAndEnergy( int numPartic
if( getNonbondedMethod() == NoCutoff || 1 ){
return calculateNoCutoffForceAndEnergy( static_cast<unsigned int>(numParticles), particlePositions, charges, dipoles, quadrupoles, tholes, dampingFactors,
polarity, axisTypes, multipoleAtomId1s, multipoleAtomId2s, forces );
polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, forces );
}
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.h
View file @ 0cae89d3
......@@ -211,8 +211,9 @@ public:
const std::vector<RealOpenMM>& dampingFactors,
const std::vector<RealOpenMM>& polarity,
const std::vector<int>& axisTypes,
const std::vector<int>& multipoleAtomId1s,
const std::vector<int>& multipoleAtomId2s,
const std::vector<int>& multipoleAtomZs,
const std::vector<int>& multipoleAtomXs,
const std::vector<int>& multipoleAtomYs,
const std::vector< std::vector< std::vector<int> > >& multipoleAtomCovalentInfo,
RealOpenMM** forces );
......@@ -461,8 +462,9 @@ private:
RealOpenMM calculateNoCutoffElectrostatic( std::vector<MultipoleParticleData>& particleData,
const std::vector<int>& axisTypes,
const std::vector<int>& multipoleAtomId1s,
const std::vector<int>& multipoleAtomId2s,
const std::vector<int>& multipoleAtomZs,
const std::vector<int>& multipoleAtomXs,
const std::vector<int>& multipoleAtomYs,
RealOpenMM** forces ) const;
/**---------------------------------------------------------------------------------------
......@@ -491,8 +493,9 @@ private:
const std::vector<RealOpenMM>& dampingFactors,
const std::vector<RealOpenMM>& polarity,
const std::vector<int>& axisTypes,
const std::vector<int>& multipoleAtomId1s,
const std::vector<int>& multipoleAtomId2s,
const std::vector<int>& multipoleAtomZs,
const std::vector<int>& multipoleAtomXs,
const std::vector<int>& multipoleAtomYs,
RealOpenMM** forces );
};
......
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