Skip to content

GitLab

Commit eb9f735a authored by leeping's avatar leeping
Browse files

Merge branch 'master' of github.com:leeping/openmm

parents 4362d539 708c4246
Show whitespace changes
Inline Side-by-side
Showing with 1562 additions and 167 deletions
.travis.yml 0 → 100644
View file @ eb9f735a
language: cpp
compiler:
- clang
before_install:
- sudo apt-get update -qq
- sudo apt-get install -qq libpcre3 libpcre3-dev gromacs
- sudo apt-get install -qq swig doxygen
- sudo apt-get install -qq python-numpy python-scipy python-nose
script:
- cmake -DCMAKE_INSTALL_PREFIX=~/OpenMM .
- make
- make test
- make install
- ls ~/OpenMM/include
- export LD_LIBRARY_PATH=~/OpenMM/lib/
- export OPENMM_LIB_PATH=~/OpenMM/lib/
- export OPENMM_INCLUDE_PATH=~/OpenMM/include/
- cd python
- sudo -E python setup.py install
- cd tests
- nosetests -vv
CMakeLists.txt
View file @ eb9f735a
......@@ -87,8 +87,11 @@ IF(WIN32)
ENDFOREACH(lib)
LINK_DIRECTORIES("${CMAKE_CURRENT_BINARY_DIR}/${CMAKE_CFG_INTDIR}")
SET(PTHREADS_LIB pthreadVC2)
SET(PTHREADS_LIB_STATIC pthreadVC2_static_mt)
ELSE(WIN32)
SET(PTHREADS_LIB pthread)
# in linux, even in static builds we link against the dynamic object (since its tied to libc versions)
SET(PTHREADS_LIB_STATIC pthread)
ENDIF(WIN32)
# The build system will set ARCH64 for 64 bit builds, which require
......
README.md
View file @ eb9f735a
## OpenMM: A High Performance Molecular Dynamics Library
[![Build Status](https://travis-ci.org/SimTk/openmm.png?branch=master)](https://travis-ci.org/SimTk/openmm)
Introduction
------------
......
openmmapi/include/OpenMM.h
View file @ eb9f735a
......@@ -62,6 +62,7 @@
#include "openmm/RBTorsionForce.h"
#include "openmm/State.h"
#include "openmm/System.h"
#include "openmm/TabulatedFunction.h"
#include "openmm/Units.h"
#include "openmm/VariableLangevinIntegrator.h"
#include "openmm/VariableVerletIntegrator.h"
......
openmmapi/include/openmm/CustomCompoundBondForce.h
View file @ eb9f735a
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -32,6 +32,7 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "TabulatedFunction.h"
#include "Force.h"
#include "Vec3.h"
#include <vector>
......@@ -91,8 +92,8 @@ namespace OpenMM {
* functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, step, delta. All trigonometric functions
* are defined in radians, and log is the natural logarithm. step(x) = 0 if x is less than 0, 1 otherwise. delta(x) = 1 if x is 0, 0 otherwise.
*
* In addition, you can call addFunction() to define a new function based on tabulated values. You specify a vector of
* values, and a natural spline is created from them. That function can then appear in the expression.
* In addition, you can call addTabulatedFunction() to define a new function based on tabulated values. You specify the function by
* creating a TabulatedFunction object. That function can then appear in the expression.
*/
class OPENMM_EXPORT CustomCompoundBondForce : public Force {
......@@ -106,6 +107,7 @@ public:
* and per-bond parameters
*/
explicit CustomCompoundBondForce(int numParticles, const std::string& energy);
~CustomCompoundBondForce();
/**
* Get the number of particles used to define each bond.
*/
......@@ -133,6 +135,14 @@ public:
/**
* Get the number of tabulated functions that have been defined.
*/
int getNumTabulatedFunctions() const {
return functions.size();
}
/**
* Get the number of tabulated functions that have been defined.
*
* @deprecated This method exists only for backward compatibility. Use getNumTabulatedFunctions() instead.
*/
int getNumFunctions() const {
return functions.size();
}
......@@ -229,33 +239,51 @@ public:
* Add a tabulated function that may appear in the energy expression.
*
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @param function a TabulatedFunction object defining the function. The TabulatedFunction
* should have been created on the heap with the "new" operator. The
* Force takes over ownership of it, and deletes it when the Force itself is deleted.
* @return the index of the function that was added
*/
int addTabulatedFunction(const std::string& name, TabulatedFunction* function);
/**
* Get a const reference to a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
const TabulatedFunction& getTabulatedFunction(int index) const;
/**
* Get a reference to a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
TabulatedFunction& getTabulatedFunction(int index);
/**
* Get the name of a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the name of the function as it appears in expressions
*/
const std::string& getTabulatedFunctionName(int index) const;
/**
* Add a tabulated function that may appear in the energy expression.
*
* @deprecated This method exists only for backward compatibility. Use addTabulatedFunction() instead.
*/
int addFunction(const std::string& name, const std::vector<double>& values, double min, double max);
/**
* Get the parameters for a tabulated function that may appear in the energy expression.
*
* @param index the index of the function for which to get parameters
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @deprecated This method exists only for backward compatibility. Use getTabulatedFunctionParameters() instead.
* If the specified function is not a Continuous1DFunction, this throws an exception.
*/
void getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const;
/**
* Set the parameters for a tabulated function that may appear in algebraic expressions.
* Set the parameters for a tabulated function that may appear in the energy expression.
*
* @param index the index of the function for which to set parameters
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @deprecated This method exists only for backward compatibility. Use setTabulatedFunctionParameters() instead.
* If the specified function is not a Continuous1DFunction, this throws an exception.
*/
void setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max);
/**
......@@ -333,12 +361,10 @@ public:
class CustomCompoundBondForce::FunctionInfo {
public:
std::string name;
std::vector<double> values;
double min, max;
TabulatedFunction* function;
FunctionInfo() {
}
FunctionInfo(const std::string& name, const std::vector<double>& values, double min, double max) :
name(name), values(values), min(min), max(max) {
FunctionInfo(const std::string& name, TabulatedFunction* function) : name(name), function(function) {
}
};
......
openmmapi/include/openmm/CustomGBForce.h
View file @ eb9f735a
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -32,6 +32,7 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "TabulatedFunction.h"
#include "Force.h"
#include "Vec3.h"
#include <map>
......@@ -134,8 +135,8 @@ namespace OpenMM {
* have the suffix "1" or "2" appended to them to indicate the values for the two interacting particles. As seen in the above example,
* an expression may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
*
* In addition, you can call addFunction() to define a new function based on tabulated values. You specify a vector of
* values, and a natural spline is created from them. That function can then appear in expressions.
* In addition, you can call addTabulatedFunction() to define a new function based on tabulated values. You specify the function by
* creating a TabulatedFunction object. That function can then appear in expressions.
*/
class OPENMM_EXPORT CustomGBForce : public Force {
......@@ -181,6 +182,7 @@ public:
* Create a CustomGBForce.
*/
CustomGBForce();
~CustomGBForce();
/**
* Get the number of particles for which force field parameters have been defined.
*/
......@@ -208,6 +210,14 @@ public:
/**
* Get the number of tabulated functions that have been defined.
*/
int getNumTabulatedFunctions() const {
return functions.size();
}
/**
* Get the number of tabulated functions that have been defined.
*
* @deprecated This method exists only for backward compatibility. Use getNumTabulatedFunctions() instead.
*/
int getNumFunctions() const {
return functions.size();
}
......@@ -452,36 +462,54 @@ public:
*/
void setExclusionParticles(int index, int particle1, int particle2);
/**
* Add a tabulated function that may appear in the energy expression.
* Add a tabulated function that may appear in expressions.
*
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @param function a TabulatedFunction object defining the function. The TabulatedFunction
* should have been created on the heap with the "new" operator. The
* Force takes over ownership of it, and deletes it when the Force itself is deleted.
* @return the index of the function that was added
*/
int addTabulatedFunction(const std::string& name, TabulatedFunction* function);
/**
* Get a const reference to a tabulated function that may appear in expressions.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
const TabulatedFunction& getTabulatedFunction(int index) const;
/**
* Get a reference to a tabulated function that may appear in expressions.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
TabulatedFunction& getTabulatedFunction(int index);
/**
* Get the name of a tabulated function that may appear in expressions.
*
* @param index the index of the function to get
* @return the name of the function as it appears in expressions
*/
const std::string& getTabulatedFunctionName(int index) const;
/**
* Add a tabulated function that may appear in expressions.
*
* @deprecated This method exists only for backward compatibility. Use addTabulatedFunction() instead.
*/
int addFunction(const std::string& name, const std::vector<double>& values, double min, double max);
/**
* Get the parameters for a tabulated function that may appear in the energy expression.
* Get the parameters for a tabulated function that may appear in expressions.
*
* @param index the index of the function for which to get parameters
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @deprecated This method exists only for backward compatibility. Use getTabulatedFunctionParameters() instead.
* If the specified function is not a Continuous1DFunction, this throws an exception.
*/
void getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const;
/**
* Set the parameters for a tabulated function that may appear in algebraic expressions.
* Set the parameters for a tabulated function that may appear in expressions.
*
* @param index the index of the function for which to set parameters
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @deprecated This method exists only for backward compatibility. Use setTabulatedFunctionParameters() instead.
* If the specified function is not a Continuous1DFunction, this throws an exception.
*/
void setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max);
/**
......@@ -577,12 +605,10 @@ public:
class CustomGBForce::FunctionInfo {
public:
std::string name;
std::vector<double> values;
double min, max;
TabulatedFunction* function;
FunctionInfo() {
}
FunctionInfo(const std::string& name, const std::vector<double>& values, double min, double max) :
name(name), values(values), min(min), max(max) {
FunctionInfo(const std::string& name, TabulatedFunction* function) : name(name), function(function) {
}
};
......
openmmapi/include/openmm/CustomHbondForce.h
View file @ eb9f735a
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -32,6 +32,7 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "TabulatedFunction.h"
#include "Force.h"
#include "Vec3.h"
#include <map>
......@@ -91,8 +92,8 @@ namespace OpenMM {
* functions: sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, tanh, erf, erfc, min, max, abs, step, delta. All trigonometric functions
* are defined in radians, and log is the natural logarithm. step(x) = 0 if x is less than 0, 1 otherwise. delta(x) = 1 if x is 0, 0 otherwise.
*
* In addition, you can call addFunction() to define a new function based on tabulated values. You specify a vector of
* values, and a natural spline is created from them. That function can then appear in the expression.
* In addition, you can call addTabulatedFunction() to define a new function based on tabulated values. You specify the function by
* creating a TabulatedFunction object. That function can then appear in the expression.
*/
class OPENMM_EXPORT CustomHbondForce : public Force {
......@@ -124,6 +125,7 @@ public:
* per-acceptor parameters
*/
explicit CustomHbondForce(const std::string& energy);
~CustomHbondForce();
/**
* Get the number of donors for which force field parameters have been defined.
*/
......@@ -163,6 +165,14 @@ public:
/**
* Get the number of tabulated functions that have been defined.
*/
int getNumTabulatedFunctions() const {
return functions.size();
}
/**
* Get the number of tabulated functions that have been defined.
*
* @deprecated This method exists only for backward compatibility. Use getNumTabulatedFunctions() instead.
*/
int getNumFunctions() const {
return functions.size();
}
......@@ -374,33 +384,51 @@ public:
* Add a tabulated function that may appear in the energy expression.
*
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @param function a TabulatedFunction object defining the function. The TabulatedFunction
* should have been created on the heap with the "new" operator. The
* Force takes over ownership of it, and deletes it when the Force itself is deleted.
* @return the index of the function that was added
*/
int addTabulatedFunction(const std::string& name, TabulatedFunction* function);
/**
* Get a const reference to a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
const TabulatedFunction& getTabulatedFunction(int index) const;
/**
* Get a reference to a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
TabulatedFunction& getTabulatedFunction(int index);
/**
* Get the name of a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the name of the function as it appears in expressions
*/
const std::string& getTabulatedFunctionName(int index) const;
/**
* Add a tabulated function that may appear in the energy expression.
*
* @deprecated This method exists only for backward compatibility. Use addTabulatedFunction() instead.
*/
int addFunction(const std::string& name, const std::vector<double>& values, double min, double max);
/**
* Get the parameters for a tabulated function that may appear in the energy expression.
*
* @param index the index of the function for which to get parameters
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @deprecated This method exists only for backward compatibility. Use getTabulatedFunctionParameters() instead.
* If the specified function is not a Continuous1DFunction, this throws an exception.
*/
void getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const;
/**
* Set the parameters for a tabulated function that may appear in algebraic expressions.
* Set the parameters for a tabulated function that may appear in the energy expression.
*
* @param index the index of the function for which to set parameters
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @deprecated This method exists only for backward compatibility. Use setTabulatedFunctionParameters() instead.
* If the specified function is not a Continuous1DFunction, this throws an exception.
*/
void setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max);
/**
......@@ -499,12 +527,10 @@ public:
class CustomHbondForce::FunctionInfo {
public:
std::string name;
std::vector<double> values;
double min, max;
TabulatedFunction* function;
FunctionInfo() {
}
FunctionInfo(const std::string& name, const std::vector<double>& values, double min, double max) :
name(name), values(values), min(min), max(max) {
FunctionInfo(const std::string& name, TabulatedFunction* function) : name(name), function(function) {
}
};
......
openmmapi/include/openmm/CustomNonbondedForce.h
View file @ eb9f735a
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2013 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -32,6 +32,7 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "TabulatedFunction.h"
#include "Force.h"
#include "Vec3.h"
#include <map>
......@@ -124,8 +125,8 @@ namespace OpenMM {
* have the suffix "1" or "2" appended to them to indicate the values for the two interacting particles. As seen in the above example,
* the expression may also involve intermediate quantities that are defined following the main expression, using ";" as a separator.
*
* In addition, you can call addFunction() to define a new function based on tabulated values. You specify a vector of
* values, and a natural spline is created from them. That function can then appear in the expression.
* In addition, you can call addTabulatedFunction() to define a new function based on tabulated values. You specify the function by
* creating a TabulatedFunction object. That function can then appear in the expression.
*/
class OPENMM_EXPORT CustomNonbondedForce : public Force {
......@@ -156,6 +157,7 @@ public:
* of r, the distance between them, as well as any global and per-particle parameters
*/
explicit CustomNonbondedForce(const std::string& energy);
~CustomNonbondedForce();
/**
* Get the number of particles for which force field parameters have been defined.
*/
......@@ -183,6 +185,14 @@ public:
/**
* Get the number of tabulated functions that have been defined.
*/
int getNumTabulatedFunctions() const {
return functions.size();
}
/**
* Get the number of tabulated functions that have been defined.
*
* @deprecated This method exists only for backward compatibility. Use getNumTabulatedFunctions() instead.
*/
int getNumFunctions() const {
return functions.size();
}
......@@ -359,33 +369,51 @@ public:
* Add a tabulated function that may appear in the energy expression.
*
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @param function a TabulatedFunction object defining the function. The TabulatedFunction
* should have been created on the heap with the "new" operator. The
* Force takes over ownership of it, and deletes it when the Force itself is deleted.
* @return the index of the function that was added
*/
int addTabulatedFunction(const std::string& name, TabulatedFunction* function);
/**
* Get a const reference to a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
const TabulatedFunction& getTabulatedFunction(int index) const;
/**
* Get a reference to a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the TabulatedFunction object defining the function
*/
TabulatedFunction& getTabulatedFunction(int index);
/**
* Get the name of a tabulated function that may appear in the energy expression.
*
* @param index the index of the function to get
* @return the name of the function as it appears in expressions
*/
const std::string& getTabulatedFunctionName(int index) const;
/**
* Add a tabulated function that may appear in the energy expression.
*
* @deprecated This method exists only for backward compatibility. Use addTabulatedFunction() instead.
*/
int addFunction(const std::string& name, const std::vector<double>& values, double min, double max);
/**
* Get the parameters for a tabulated function that may appear in the energy expression.
*
* @param index the index of the function for which to get parameters
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @deprecated This method exists only for backward compatibility. Use getTabulatedFunctionParameters() instead.
* If the specified function is not a Continuous1DFunction, this throws an exception.
*/
void getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const;
/**
* Set the parameters for a tabulated function that may appear in algebraic expressions.
* Set the parameters for a tabulated function that may appear in the energy expression.
*
* @param index the index of the function for which to set parameters
* @param name the name of the function as it appears in expressions
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min and max.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of the independent variable corresponding to the first element of values
* @param max the value of the independent variable corresponding to the last element of values
* @deprecated This method exists only for backward compatibility. Use setTabulatedFunctionParameters() instead.
* If the specified function is not a Continuous1DFunction, this throws an exception.
*/
void setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max);
/**
......@@ -507,12 +535,10 @@ public:
class CustomNonbondedForce::FunctionInfo {
public:
std::string name;
std::vector<double> values;
double min, max;
TabulatedFunction* function;
FunctionInfo() {
}
FunctionInfo(const std::string& name, const std::vector<double>& values, double min, double max) :
name(name), values(values), min(min), max(max) {
FunctionInfo(const std::string& name, TabulatedFunction* function) : name(name), function(function) {
}
};
......
openmmapi/include/openmm/TabulatedFunction.h 0 → 100644
View file @ eb9f735a
#ifndef OPENMM_TABULATEDFUNCTION_H_
#define OPENMM_TABULATEDFUNCTION_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "internal/windowsExport.h"
#include <vector>
namespace OpenMM {
/**
* A TabulatedFunction uses a set of tabulated values to define a mathematical function.
* It can be used by various custom forces.
*
* TabulatedFunction is an abstract class with concrete subclasses for more specific
* types of functions. There are subclasses for:
*
* <ul>
* <li>1, 2, and 3 dimensional functions. The dimensionality of a function means
* the number of input arguments it takes.</li>
* <li>Continuous and discrete functions. A continuous function is interpolated by
* fitting a natural cubic spline to the tabulated values. A discrete function is
* only defined for integer values of its arguments (that is, at the tabulated points),
* and does not try to interpolate between them. Discrete function can be evaluated
* more quickly than continuous ones.</li>
* </ul>
*/
class OPENMM_EXPORT TabulatedFunction {
public:
virtual ~TabulatedFunction() {
}
};
/**
* This is a TabulatedFunction that computes a continuous one dimensional function.
*/
class OPENMM_EXPORT Continuous1DFunction : public TabulatedFunction {
public:
/**
* Create a Continuous1DFunction f(x) based on a set of tabulated values.
*
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min
* and max. A natural cubic spline is used to interpolate between the tabulated values.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of x corresponding to the first element of values
* @param max the value of x corresponding to the last element of values
*/
Continuous1DFunction(const std::vector<double>& values, double min, double max);
/**
* Get the parameters for the tabulated function.
*
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min
* and max. A natural cubic spline is used to interpolate between the tabulated values.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of x corresponding to the first element of values
* @param max the value of x corresponding to the last element of values
*/
void getFunctionParameters(std::vector<double>& values, double& min, double& max) const;
/**
* Set the parameters for the tabulated function.
*
* @param values the tabulated values of the function f(x) at uniformly spaced values of x between min
* and max. A natural cubic spline is used to interpolate between the tabulated values.
* The function is assumed to be zero for x &lt; min or x &gt; max.
* @param min the value of x corresponding to the first element of values
* @param max the value of x corresponding to the last element of values
*/
void setFunctionParameters(const std::vector<double>& values, double min, double max);
private:
std::vector<double> values;
double min, max;
};
/**
* This is a TabulatedFunction that computes a continuous two dimensional function.
*/
class OPENMM_EXPORT Continuous2DFunction : public TabulatedFunction {
public:
/**
* Create a Continuous2DFunction f(x,y) based on a set of tabulated values.
*
* @param values the tabulated values of the function f(x,y) at xsize uniformly spaced values of x between xmin
* and xmax, and ysize values of y between ymin and ymax. A natural cubic spline is used to interpolate between the tabulated values.
* The function is assumed to be zero when x or y is outside its specified range. The values should be ordered so that
* values[i+xsize*j] = f(x_i,y_j), where x_i is the i'th uniformly spaced value of x. This must be of length xsize*ysize.
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param xmin the value of x corresponding to the first element of values
* @param xmax the value of x corresponding to the last element of values
* @param ymin the value of y corresponding to the first element of values
* @param ymax the value of y corresponding to the last element of values
*/
Continuous2DFunction(int xsize, int ysize, const std::vector<double>& values, double xmin, double xmax, double ymin, double ymax);
/**
* Get the parameters for the tabulated function.
*
* @param values the tabulated values of the function f(x,y) at xsize uniformly spaced values of x between xmin
* and xmax, and ysize values of y between ymin and ymax. A natural cubic spline is used to interpolate between the tabulated values.
* The function is assumed to be zero when x or y is outside its specified range. The values should be ordered so that
* values[i+xsize*j] = f(x_i,y_j), where x_i is the i'th uniformly spaced value of x. This must be of length xsize*ysize.
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param xmin the value of x corresponding to the first element of values
* @param xmax the value of x corresponding to the last element of values
* @param ymin the value of y corresponding to the first element of values
* @param ymax the value of y corresponding to the last element of values
*/
void getFunctionParameters(int& xsize, int& ysize, std::vector<double>& values, double& xmin, double& xmax, double& ymin, double& ymax) const;
/**
* Set the parameters for the tabulated function.
*
* @param values the tabulated values of the function f(x,y) at xsize uniformly spaced values of x between xmin
* and xmax, and ysize values of y between ymin and ymax. A natural cubic spline is used to interpolate between the tabulated values.
* The function is assumed to be zero when x or y is outside its specified range. The values should be ordered so that
* values[i+xsize*j] = f(x_i,y_j), where x_i is the i'th uniformly spaced value of x. This must be of length xsize*ysize.
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param xmin the value of x corresponding to the first element of values
* @param xmax the value of x corresponding to the last element of values
* @param ymin the value of y corresponding to the first element of values
* @param ymax the value of y corresponding to the last element of values
*/
void setFunctionParameters(int xsize, int ysize, const std::vector<double>& values, double xmin, double xmax, double ymin, double ymax);
private:
std::vector<double> values;
int xsize, ysize;
double xmin, xmax, ymin, ymax;
};
/**
* This is a TabulatedFunction that computes a continuous three dimensional function.
*/
class OPENMM_EXPORT Continuous3DFunction : public TabulatedFunction {
public:
/**
* Create a Continuous3DFunction f(x,y,z) based on a set of tabulated values.
*
* @param values the tabulated values of the function f(x,y,z) at xsize uniformly spaced values of x between xmin
* and xmax, ysize values of y between ymin and ymax, and zsize values of z between zmin and zmax.
* A natural cubic spline is used to interpolate between the tabulated values. The function is
* assumed to be zero when x, y, or z is outside its specified range. The values should be ordered so
* that values[i+xsize*j+xsize*ysize*k] = f(x_i,y_j,z_k), where x_i is the i'th uniformly spaced value of x.
* This must be of length xsize*ysize*zsize.
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param ysize the number of table elements along the z direction
* @param xmin the value of x corresponding to the first element of values
* @param xmax the value of x corresponding to the last element of values
* @param ymin the value of y corresponding to the first element of values
* @param ymax the value of y corresponding to the last element of values
* @param zmin the value of z corresponding to the first element of values
* @param zmax the value of z corresponding to the last element of values
*/
Continuous3DFunction(int xsize, int ysize, int zsize, const std::vector<double>& values, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax);
/**
* Get the parameters for the tabulated function.
*
* @param values the tabulated values of the function f(x,y,z) at xsize uniformly spaced values of x between xmin
* and xmax, ysize values of y between ymin and ymax, and zsize values of z between zmin and zmax.
* A natural cubic spline is used to interpolate between the tabulated values. The function is
* assumed to be zero when x, y, or z is outside its specified range. The values should be ordered so
* that values[i+xsize*j+xsize*ysize*k] = f(x_i,y_j,z_k), where x_i is the i'th uniformly spaced value of x.
* This must be of length xsize*ysize*zsize.
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param ysize the number of table elements along the z direction
* @param xmin the value of x corresponding to the first element of values
* @param xmax the value of x corresponding to the last element of values
* @param ymin the value of y corresponding to the first element of values
* @param ymax the value of y corresponding to the last element of values
* @param zmin the value of z corresponding to the first element of values
* @param zmax the value of z corresponding to the last element of values
*/
void getFunctionParameters(int& xsize, int& ysize, int& zsize, std::vector<double>& values, double& xmin, double& xmax, double& ymin, double& ymax, double& zmin, double& zmax) const;
/**
* Set the parameters for the tabulated function.
*
* @param values the tabulated values of the function f(x,y,z) at xsize uniformly spaced values of x between xmin
* and xmax, ysize values of y between ymin and ymax, and zsize values of z between zmin and zmax.
* A natural cubic spline is used to interpolate between the tabulated values. The function is
* assumed to be zero when x, y, or z is outside its specified range. The values should be ordered so
* that values[i+xsize*j+xsize*ysize*k] = f(x_i,y_j,z_k), where x_i is the i'th uniformly spaced value of x.
* This must be of length xsize*ysize*zsize.
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param ysize the number of table elements along the z direction
* @param xmin the value of x corresponding to the first element of values
* @param xmax the value of x corresponding to the last element of values
* @param ymin the value of y corresponding to the first element of values
* @param ymax the value of y corresponding to the last element of values
* @param zmin the value of z corresponding to the first element of values
* @param zmax the value of z corresponding to the last element of values
*/
void setFunctionParameters(int xsize, int ysize, int zsize, const std::vector<double>& values, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax);
private:
std::vector<double> values;
int xsize, ysize, zsize;
double xmin, xmax, ymin, ymax, zmin, zmax;
};
/**
* This is a TabulatedFunction that computes a discrete one dimensional function f(x).
* To evaluate it, x is rounded to the nearest integer and the table element with that
* index is returned. If the index is outside the range [0, size), the result is undefined.
*/
class OPENMM_EXPORT Discrete1DFunction : public TabulatedFunction {
public:
/**
* Create a Discrete1DFunction f(x) based on a set of tabulated values.
*
* @param values the tabulated values of the function f(x)
*/
Discrete1DFunction(const std::vector<double>& values);
/**
* Get the parameters for the tabulated function.
*
* @param values the tabulated values of the function f(x)
*/
void getFunctionParameters(std::vector<double>& values) const;
/**
* Set the parameters for the tabulated function.
*
* @param values the tabulated values of the function f(x)
*/
void setFunctionParameters(const std::vector<double>& values);
private:
std::vector<double> values;
};
/**
* This is a TabulatedFunction that computes a discrete two dimensional function f(x,y).
* To evaluate it, x and y are each rounded to the nearest integer and the table element with those
* indices is returned. If either index is outside the range [0, size), the result is undefined.
*/
class OPENMM_EXPORT Discrete2DFunction : public TabulatedFunction {
public:
/**
* Create a Discrete2DFunction f(x,y) based on a set of tabulated values.
*
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param values the tabulated values of the function f(x,y), ordered so that
* values[i+xsize*j] = f(i,j). This must be of length xsize*ysize.
*/
Discrete2DFunction(int xsize, int ysize, const std::vector<double>& values);
/**
* Get the parameters for the tabulated function.
*
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param values the tabulated values of the function f(x,y), ordered so that
* values[i+xsize*j] = f(i,j). This must be of length xsize*ysize.
*/
void getFunctionParameters(int& xsize, int& ysize, std::vector<double>& values) const;
/**
* Set the parameters for the tabulated function.
*
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param values the tabulated values of the function f(x,y), ordered so that
* values[i+xsize*j] = f(i,j). This must be of length xsize*ysize.
*/
void setFunctionParameters(int xsize, int ysize, const std::vector<double>& values);
private:
int xsize, ysize;
std::vector<double> values;
};
/**
* This is a TabulatedFunction that computes a discrete three dimensional function f(x,y,z).
* To evaluate it, x, y, and z are each rounded to the nearest integer and the table element with those
* indices is returned. If any index is outside the range [0, size), the result is undefined.
*/
class OPENMM_EXPORT Discrete3DFunction : public TabulatedFunction {
public:
/**
* Create a Discrete3DFunction f(x,y,z) based on a set of tabulated values.
*
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param zsize the number of table elements along the z direction
* @param values the tabulated values of the function f(x,y,z), ordered so that
* values[i+xsize*j+xsize*ysize*k] = f(i,j,k). This must be of length xsize*ysize*zsize.
*/
Discrete3DFunction(int xsize, int ysize, int zsize, const std::vector<double>& values);
/**
* Get the parameters for the tabulated function.
*
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param zsize the number of table elements along the z direction
* @param values the tabulated values of the function f(x,y,z), ordered so that
* values[i+xsize*j+xsize*ysize*k] = f(i,j,k). This must be of length xsize*ysize*zsize.
*/
void getFunctionParameters(int& xsize, int& ysize, int& zsize, std::vector<double>& values) const;
/**
* Set the parameters for the tabulated function.
*
* @param xsize the number of table elements along the x direction
* @param ysize the number of table elements along the y direction
* @param zsize the number of table elements along the z direction
* @param values the tabulated values of the function f(x,y,z), ordered so that
* values[i+xsize*j+xsize*ysize*k] = f(i,j,k). This must be of length xsize*ysize*zsize.
*/
void setFunctionParameters(int xsize, int ysize, int zsize, const std::vector<double>& values);
private:
int xsize, ysize, zsize;
std::vector<double> values;
};
} // namespace OpenMM
#endif /*OPENMM_TABULATEDFUNCTION_H_*/
openmmapi/include/openmm/internal/SplineFitter.h
View file @ eb9f735a
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2010 Stanford University and the Authors. *
* Portions copyright (c) 2010-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -67,7 +67,7 @@ public:
*/
static void createPeriodicSpline(const std::vector<double>& x, const std::vector<double>& y, std::vector<double>& deriv);
/**
* Evaluate a spline generated by one of the other methods in this class.
* Evaluate a 1D spline generated by one of the other methods in this class.
*
* @param x the values of the independent variable at the data points to interpolate
* @param y the values of the dependent variable at the data points to interpolate
......@@ -77,7 +77,7 @@ public:
*/
static double evaluateSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& deriv, double t);
/**
* Evaluate the derivative of a spline generated by one of the other methods in this class.
* Evaluate the derivative of a 1D spline generated by one of the other methods in this class.
*
* @param x the values of the independent variable at the data points to interpolate
* @param y the values of the dependent variable at the data points to interpolate
......@@ -86,6 +86,90 @@ public:
* @return the value of the spline's derivative at the specified point
*/
static double evaluateSplineDerivative(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& deriv, double t);
/**
* Fit a natural cubic spline surface f(x,y) to a 2D set of data points. The resulting spline interpolates all the
* data points, has a continuous second derivative everywhere, and has a second derivative of 0 at the boundary.
*
* @param x the values of the first independent variable at the data points to interpolate. They must
* be strictly increasing: x[i] > x[i-1].
* @param y the values of the second independent variable at the data points to interpolate. They must
* be strictly increasing: y[i] > y[i-1].
* @param values the values of the dependent variable at the data points to interpolate. They must be ordered
* so that values[i+xsize*j] = f(x[i],y[j]), where xsize is the length of x.
* @param c on exit, this contains the spline coefficients at each of the data points
*/
static void create2DNaturalSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& values, std::vector<std::vector<double> >& c);
/**
* Evaluate a 2D spline generated by one of the other methods in this class.
*
* @param x the values of the first independent variable at the data points to interpolate
* @param y the values of the second independent variable at the data points to interpolate
* @param values the values of the dependent variable at the data points to interpolate
* @param c the vector of spline coefficients that was calculated by one of the other methods
* @param u the value of the first independent variable at which to evaluate the spline
* @param v the value of the second independent variable at which to evaluate the spline
* @return the value of the spline at the specified point
*/
static double evaluate2DSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& values, const std::vector<std::vector<double> >& c, double u, double v);
/**
* Evaluate the derivatives of a 2D spline generated by one of the other methods in this class.
*
* @param x the values of the first independent variable at the data points to interpolate
* @param y the values of the second independent variable at the data points to interpolate
* @param values the values of the dependent variable at the data points to interpolate
* @param c the vector of spline coefficients that was calculated by one of the other methods
* @param u the value of the first independent variable at which to evaluate the spline
* @param v the value of the second independent variable at which to evaluate the spline
* @param dx on exit, the x derivative of the spline at the specified point
* @param dy on exit, the y derivative of the spline at the specified point
*/
static void evaluate2DSplineDerivatives(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& values, const std::vector<std::vector<double> >& c, double u, double v, double& dx, double& dy);
/**
* Fit a natural cubic spline surface f(x,y,z) to a 3D set of data points. The resulting spline interpolates all the
* data points, has a continuous second derivative everywhere, and has a second derivative of 0 at the boundary.
*
* @param x the values of the first independent variable at the data points to interpolate. They must
* be strictly increasing: x[i] > x[i-1].
* @param y the values of the second independent variable at the data points to interpolate. They must
* be strictly increasing: y[i] > y[i-1].
* @param z the values of the third independent variable at the data points to interpolate. They must
* be strictly increasing: z[i] > z[i-1].
* @param values the values of the dependent variable at the data points to interpolate. They must be ordered
* so that values[i+xsize*j+xsize*ysize*k] = f(x[i],y[j],z[k]), where xsize is the length of x
* and ysize is the length of y.
* @param c on exit, this contains the spline coefficients at each of the data points
*/
static void create3DNaturalSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, const std::vector<double>& values, std::vector<std::vector<double> >& c);
/**
* Evaluate a 3D spline generated by one of the other methods in this class.
*
* @param x the values of the first independent variable at the data points to interpolate
* @param y the values of the second independent variable at the data points to interpolate
* @param z the values of the third independent variable at the data points to interpolate
* @param values the values of the dependent variable at the data points to interpolate
* @param c the vector of spline coefficients that was calculated by one of the other methods
* @param u the value of the first independent variable at which to evaluate the spline
* @param v the value of the second independent variable at which to evaluate the spline
* @param w the value of the third independent variable at which to evaluate the spline
* @return the value of the spline at the specified point
*/
static double evaluate3DSpline(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, const std::vector<double>& values, const std::vector<std::vector<double> >& c, double u, double v, double w);
/**
* Evaluate the derivatives of a 3D spline generated by one of the other methods in this class.
*
* @param x the values of the first independent variable at the data points to interpolate
* @param y the values of the second independent variable at the data points to interpolate
* @param z the values of the third independent variable at the data points to interpolate
* @param values the values of the dependent variable at the data points to interpolate
* @param c the vector of spline coefficients that was calculated by one of the other methods
* @param u the value of the first independent variable at which to evaluate the spline
* @param v the value of the second independent variable at which to evaluate the spline
* @param w the value of the third independent variable at which to evaluate the spline
* @param dx on exit, the x derivative of the spline at the specified point
* @param dy on exit, the y derivative of the spline at the specified point
* @param dz on exit, the z derivative of the spline at the specified point
*/
static void evaluate3DSplineDerivatives(const std::vector<double>& x, const std::vector<double>& y, const std::vector<double>& z, const std::vector<double>& values, const std::vector<std::vector<double> >& c, double u, double v, double w, double& dx, double& dy, double &dz);
private:
static void solveTridiagonalMatrix(const std::vector<double>& a, const std::vector<double>& b, const std::vector<double>& c, const std::vector<double>& rhs, std::vector<double>& sol);
};
......
openmmapi/src/CustomCompoundBondForce.cpp
View file @ eb9f735a
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -51,6 +51,12 @@ using std::vector;
CustomCompoundBondForce::CustomCompoundBondForce(int numParticles, const string& energy) : particlesPerBond(numParticles), energyExpression(energy) {
}
CustomCompoundBondForce::~CustomCompoundBondForce() {
for (int i = 0; i < (int) functions.size(); i++)
delete functions[i].function;
}
const string& CustomCompoundBondForce::getEnergyFunction() const {
return energyExpression;
}
......@@ -120,33 +126,47 @@ void CustomCompoundBondForce::setBondParameters(int index, const vector<int>& pa
bonds[index].parameters = parameters;
}
int CustomCompoundBondForce::addTabulatedFunction(const std::string& name, TabulatedFunction* function) {
functions.push_back(FunctionInfo(name, function));
return functions.size()-1;
}
const TabulatedFunction& CustomCompoundBondForce::getTabulatedFunction(int index) const {
ASSERT_VALID_INDEX(index, functions);
return *functions[index].function;
}
TabulatedFunction& CustomCompoundBondForce::getTabulatedFunction(int index) {
ASSERT_VALID_INDEX(index, functions);
return *functions[index].function;
}
const string& CustomCompoundBondForce::getTabulatedFunctionName(int index) const {
ASSERT_VALID_INDEX(index, functions);
return functions[index].name;
}
int CustomCompoundBondForce::addFunction(const std::string& name, const std::vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("CustomCompoundBondForce: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("CustomCompoundBondForce: a tabulated function must have at least two points");
functions.push_back(FunctionInfo(name, values, min, max));
functions.push_back(FunctionInfo(name, new Continuous1DFunction(values, min, max)));
return functions.size()-1;
}
void CustomCompoundBondForce::getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const {
ASSERT_VALID_INDEX(index, functions);
Continuous1DFunction* function = dynamic_cast<Continuous1DFunction*>(functions[index].function);
if (function == NULL)
throw OpenMMException("CustomCompoundBondForce: function is not a Continuous1DFunction");
name = functions[index].name;
values = functions[index].values;
min = functions[index].min;
max = functions[index].max;
function->getFunctionParameters(values, min, max);
}
void CustomCompoundBondForce::setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("CustomCompoundBondForce: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("CustomCompoundBondForce: a tabulated function must have at least two points");
ASSERT_VALID_INDEX(index, functions);
Continuous1DFunction* function = dynamic_cast<Continuous1DFunction*>(functions[index].function);
if (function == NULL)
throw OpenMMException("CustomCompoundBondForce: function is not a Continuous1DFunction");
functions[index].name = name;
functions[index].values = values;
functions[index].min = min;
functions[index].max = max;
function->setFunctionParameters(values, min, max);
}
ForceImpl* CustomCompoundBondForce::createImpl() const {
......
openmmapi/src/CustomCompoundBondForceImpl.cpp
View file @ eb9f735a
......@@ -147,7 +147,7 @@ ParsedExpression CustomCompoundBondForceImpl::prepareExpression(const CustomComp
ExpressionTreeNode CustomCompoundBondForceImpl::replaceFunctions(const ExpressionTreeNode& node, map<string, int> atoms,
map<string, vector<int> >& distances, map<string, vector<int> >& angles, map<string, vector<int> >& dihedrals) {
const Operation& op = node.getOperation();
if (op.getId() != Operation::CUSTOM || op.getNumArguments() < 2)
if (op.getId() != Operation::CUSTOM || (op.getName() != "distance" && op.getName() != "angle" && op.getName() != "dihedral"))
{
// This is not an angle or dihedral, so process its children.
......
openmmapi/src/CustomGBForce.cpp
View file @ eb9f735a
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -50,6 +50,11 @@ using std::vector;
CustomGBForce::CustomGBForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0) {
}
CustomGBForce::~CustomGBForce() {
for (int i = 0; i < (int) functions.size(); i++)
delete functions[i].function;
}
CustomGBForce::NonbondedMethod CustomGBForce::getNonbondedMethod() const {
return nonbondedMethod;
}
......@@ -173,33 +178,47 @@ void CustomGBForce::setExclusionParticles(int index, int particle1, int particle
exclusions[index].particle2 = particle2;
}
int CustomGBForce::addTabulatedFunction(const std::string& name, TabulatedFunction* function) {
functions.push_back(FunctionInfo(name, function));
return functions.size()-1;
}
const TabulatedFunction& CustomGBForce::getTabulatedFunction(int index) const {
ASSERT_VALID_INDEX(index, functions);
return *functions[index].function;
}
TabulatedFunction& CustomGBForce::getTabulatedFunction(int index) {
ASSERT_VALID_INDEX(index, functions);
return *functions[index].function;
}
const string& CustomGBForce::getTabulatedFunctionName(int index) const {
ASSERT_VALID_INDEX(index, functions);
return functions[index].name;
}
int CustomGBForce::addFunction(const std::string& name, const std::vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("CustomGBForce: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("CustomGBForce: a tabulated function must have at least two points");
functions.push_back(FunctionInfo(name, values, min, max));
functions.push_back(FunctionInfo(name, new Continuous1DFunction(values, min, max)));
return functions.size()-1;
}
void CustomGBForce::getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const {
ASSERT_VALID_INDEX(index, functions);
Continuous1DFunction* function = dynamic_cast<Continuous1DFunction*>(functions[index].function);
if (function == NULL)
throw OpenMMException("CustomGBForce: function is not a Continuous1DFunction");
name = functions[index].name;
values = functions[index].values;
min = functions[index].min;
max = functions[index].max;
function->getFunctionParameters(values, min, max);
}
void CustomGBForce::setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("CustomGBForce: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("CustomGBForce: a tabulated function must have at least two points");
ASSERT_VALID_INDEX(index, functions);
Continuous1DFunction* function = dynamic_cast<Continuous1DFunction*>(functions[index].function);
if (function == NULL)
throw OpenMMException("CustomGBForce: function is not a Continuous1DFunction");
functions[index].name = name;
functions[index].values = values;
functions[index].min = min;
functions[index].max = max;
function->setFunctionParameters(values, min, max);
}
ForceImpl* CustomGBForce::createImpl() const {
......
openmmapi/src/CustomHbondForce.cpp
View file @ eb9f735a
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -50,6 +50,12 @@ using std::vector;
CustomHbondForce::CustomHbondForce(const string& energy) : energyExpression(energy), nonbondedMethod(NoCutoff), cutoffDistance(1.0) {
}
CustomHbondForce::~CustomHbondForce() {
for (int i = 0; i < (int) functions.size(); i++)
delete functions[i].function;
}
const string& CustomHbondForce::getEnergyFunction() const {
return energyExpression;
}
......@@ -187,33 +193,47 @@ void CustomHbondForce::setExclusionParticles(int index, int donor, int acceptor)
exclusions[index].acceptor = acceptor;
}
int CustomHbondForce::addTabulatedFunction(const std::string& name, TabulatedFunction* function) {
functions.push_back(FunctionInfo(name, function));
return functions.size()-1;
}
const TabulatedFunction& CustomHbondForce::getTabulatedFunction(int index) const {
ASSERT_VALID_INDEX(index, functions);
return *functions[index].function;
}
TabulatedFunction& CustomHbondForce::getTabulatedFunction(int index) {
ASSERT_VALID_INDEX(index, functions);
return *functions[index].function;
}
const string& CustomHbondForce::getTabulatedFunctionName(int index) const {
ASSERT_VALID_INDEX(index, functions);
return functions[index].name;
}
int CustomHbondForce::addFunction(const std::string& name, const std::vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("CustomHbondForce: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("CustomHbondForce: a tabulated function must have at least two points");
functions.push_back(FunctionInfo(name, values, min, max));
functions.push_back(FunctionInfo(name, new Continuous1DFunction(values, min, max)));
return functions.size()-1;
}
void CustomHbondForce::getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const {
ASSERT_VALID_INDEX(index, functions);
Continuous1DFunction* function = dynamic_cast<Continuous1DFunction*>(functions[index].function);
if (function == NULL)
throw OpenMMException("CustomHbondForce: function is not a Continuous1DFunction");
name = functions[index].name;
values = functions[index].values;
min = functions[index].min;
max = functions[index].max;
function->getFunctionParameters(values, min, max);
}
void CustomHbondForce::setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("CustomHbondForce: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("CustomHbondForce: a tabulated function must have at least two points");
ASSERT_VALID_INDEX(index, functions);
Continuous1DFunction* function = dynamic_cast<Continuous1DFunction*>(functions[index].function);
if (function == NULL)
throw OpenMMException("CustomHbondForce: function is not a Continuous1DFunction");
functions[index].name = name;
functions[index].values = values;
functions[index].min = min;
functions[index].max = max;
function->setFunctionParameters(values, min, max);
}
ForceImpl* CustomHbondForce::createImpl() const {
......
openmmapi/src/CustomNonbondedForce.cpp
View file @ eb9f735a
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -51,6 +51,11 @@ CustomNonbondedForce::CustomNonbondedForce(const string& energy) : energyExpress
switchingDistance(-1.0), useSwitchingFunction(false), useLongRangeCorrection(false) {
}
CustomNonbondedForce::~CustomNonbondedForce() {
for (int i = 0; i < (int) functions.size(); i++)
delete functions[i].function;
}
const string& CustomNonbondedForce::getEnergyFunction() const {
return energyExpression;
}
......@@ -169,34 +174,47 @@ void CustomNonbondedForce::setExclusionParticles(int index, int particle1, int p
exclusions[index].particle1 = particle1;
exclusions[index].particle2 = particle2;
}
int CustomNonbondedForce::addTabulatedFunction(const std::string& name, TabulatedFunction* function) {
functions.push_back(FunctionInfo(name, function));
return functions.size()-1;
}
const TabulatedFunction& CustomNonbondedForce::getTabulatedFunction(int index) const {
ASSERT_VALID_INDEX(index, functions);
return *functions[index].function;
}
TabulatedFunction& CustomNonbondedForce::getTabulatedFunction(int index) {
ASSERT_VALID_INDEX(index, functions);
return *functions[index].function;
}
const string& CustomNonbondedForce::getTabulatedFunctionName(int index) const {
ASSERT_VALID_INDEX(index, functions);
return functions[index].name;
}
int CustomNonbondedForce::addFunction(const std::string& name, const std::vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("CustomNonbondedForce: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("CustomNonbondedForce: a tabulated function must have at least two points");
functions.push_back(FunctionInfo(name, values, min, max));
functions.push_back(FunctionInfo(name, new Continuous1DFunction(values, min, max)));
return functions.size()-1;
}
void CustomNonbondedForce::getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const {
ASSERT_VALID_INDEX(index, functions);
Continuous1DFunction* function = dynamic_cast<Continuous1DFunction*>(functions[index].function);
if (function == NULL)
throw OpenMMException("CustomNonbondedForce: function is not a Continuous1DFunction");
name = functions[index].name;
values = functions[index].values;
min = functions[index].min;
max = functions[index].max;
function->getFunctionParameters(values, min, max);
}
void CustomNonbondedForce::setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("CustomNonbondedForce: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("CustomNonbondedForce: a tabulated function must have at least two points");
ASSERT_VALID_INDEX(index, functions);
Continuous1DFunction* function = dynamic_cast<Continuous1DFunction*>(functions[index].function);
if (function == NULL)
throw OpenMMException("CustomNonbondedForce: function is not a Continuous1DFunction");
functions[index].name = name;
functions[index].values = values;
functions[index].min = min;
functions[index].max = max;
function->setFunctionParameters(values, min, max);
}
int CustomNonbondedForce::addInteractionGroup(const std::set<int>& set1, const std::set<int>& set2) {
......
openmmapi/src/SplineFitter.cpp
View file @ eb9f735a
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2010 Stanford University and the Authors. *
* Portions copyright (c) 2010-2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -190,3 +190,535 @@ void SplineFitter::solveTridiagonalMatrix(const vector<double>& a, const vector<
for (int i = n-2; i >= 0; i--)
sol[i] -= gamma[i+1]*sol[i+1];
}
void SplineFitter::create2DNaturalSpline(const vector<double>& x, const vector<double>& y, const vector<double>& values, vector<vector<double> >& c) {
int xsize = x.size(), ysize = y.size();
if (xsize < 2 || ysize < 2)
throw OpenMMException("create2DNaturalSpline: must have at least two points along each axis");
if (values.size() != xsize*ysize)
throw OpenMMException("create2DNaturalSpline: incorrect number of values");
vector<double> d1(xsize*ysize), d2(xsize*ysize), d12(xsize*ysize);
vector<double> t(xsize), deriv(xsize);
// Compute derivatives with respect to x.
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < xsize; j++)
t[j] = values[j+xsize*i];
SplineFitter::createNaturalSpline(x, t, deriv);
for (int j = 0; j < xsize; j++)
d1[j+xsize*i] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[j]);
}
// Compute derivatives with respect to y.
t.resize(ysize);
deriv.resize(ysize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++)
t[j] = values[i+xsize*j];
SplineFitter::createNaturalSpline(y, t, deriv);
for (int j = 0; j < ysize; j++)
d2[i+xsize*j] = SplineFitter::evaluateSplineDerivative(y, t, deriv, y[j]);
}
// Compute cross derivatives.
t.resize(xsize);
deriv.resize(xsize);
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < xsize; j++)
t[j] = d2[j+xsize*i];
SplineFitter::createNaturalSpline(x, t, deriv);
for (int j = 0; j < xsize; j++)
d12[j+xsize*i] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[j]);
}
// Now compute the coefficients.
const int wt[] = {
1, 0, -3, 2, 0, 0, 0, 0, -3, 0, 9, -6, 2, 0, -6, 4,
0, 0, 0, 0, 0, 0, 0, 0, 3, 0, -9, 6, -2, 0, 6, -4,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -6, 0, 0, -6, 4,
0, 0, 3, -2, 0, 0, 0, 0, 0, 0, -9, 6, 0, 0, 6, -4,
0, 0, 0, 0, 1, 0, -3, 2, -2, 0, 6, -4, 1, 0, -3, 2,
0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 3, -2, 1, 0, -3, 2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 2, 0, 0, 3, -2,
0, 0, 0, 0, 0, 0, 3, -2, 0, 0, -6, 4, 0, 0, 3, -2,
0, 1, -2, 1, 0, 0, 0, 0, 0, -3, 6, -3, 0, 2, -4, 2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 3, -6, 3, 0, -2, 4, -2,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 2, -2,
0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 3, -3, 0, 0, -2, 2,
0, 0, 0, 0, 0, 1, -2, 1, 0, -2, 4, -2, 0, 1, -2, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 2, -1, 0, 1, -2, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0, 0, -1, 1,
0, 0, 0, 0, 0, 0, -1, 1, 0, 0, 2, -2, 0, 0, -1, 1
};
vector<double> rhs(16);
c.resize((xsize-1)*(ysize-1));
for (int i = 0; i < xsize-1; i++) {
for (int j = 0; j < ysize-1; j++) {
// Compute the 16 coefficients for patch (i, j).
int nexti = i+1;
int nextj = j+1;
double deltax = x[nexti]-x[i];
double deltay = y[nextj]-y[j];
double e[] = {values[i+j*xsize], values[nexti+j*xsize], values[nexti+nextj*xsize], values[i+nextj*xsize]};
double e1[] = {d1[i+j*xsize], d1[nexti+j*xsize], d1[nexti+nextj*xsize], d1[i+nextj*xsize]};
double e2[] = {d2[i+j*xsize], d2[nexti+j*xsize], d2[nexti+nextj*xsize], d2[i+nextj*xsize]};
double e12[] = {d12[i+j*xsize], d12[nexti+j*xsize], d12[nexti+nextj*xsize], d12[i+nextj*xsize]};
for (int k = 0; k < 4; k++) {
rhs[k] = e[k];
rhs[k+4] = e1[k]*deltax;
rhs[k+8] = e2[k]*deltay;
rhs[k+12] = e12[k]*deltax*deltay;
}
vector<double>& coeff = c[i+j*(xsize-1)];
coeff.resize(16);
for (int k = 0; k < 16; k++) {
double sum = 0.0;
for (int m = 0; m < 16; m++)
sum += wt[k+16*m]*rhs[m];
coeff[k] = sum;
}
}
}
}
double SplineFitter::evaluate2DSpline(const vector<double>& x, const vector<double>& y, const vector<double>& values, const vector<vector<double> >& c, double u, double v) {
int xsize = x.size();
int ysize = y.size();
if (u < x[0] || u > x[xsize-1] || v < y[0] || v > y[ysize-1])
throw OpenMMException("evaluate2DSpline: specified point is outside the range defined by the spline");
// Perform a binary search to identify the interval containing the point to evaluate.
int lowerx = 0;
int upperx = xsize-1;
while (upperx-lowerx > 1) {
int middle = (upperx+lowerx)/2;
if (x[middle] > u)
upperx = middle;
else
lowerx = middle;
}
int lowery = 0;
int uppery = ysize-1;
while (uppery-lowery > 1) {
int middle = (uppery+lowery)/2;
if (y[middle] > v)
uppery = middle;
else
lowery = middle;
}
double deltax = x[upperx]-x[lowerx];
double deltay = y[uppery]-y[lowery];
double da = (u-x[lowerx])/deltax;
double db = (v-y[lowery])/deltay;
const vector<double>& coeff = c[lowerx+(xsize-1)*lowery];
// Evaluate the spline to determine the value.
double value = 0;
for (int i = 3; i >= 0; i--)
value = da*value + ((coeff[i*4+3]*db + coeff[i*4+2])*db + coeff[i*4+1])*db + coeff[i*4+0];
return value;
}
void SplineFitter::evaluate2DSplineDerivatives(const vector<double>& x, const vector<double>& y, const vector<double>& values, const vector<vector<double> >& c, double u, double v, double& dx, double &dy) {
int xsize = x.size();
int ysize = y.size();
if (u < x[0] || u > x[xsize-1] || v < y[0] || v > y[ysize-1])
throw OpenMMException("evaluate2DSplineDerivatives: specified point is outside the range defined by the spline");
// Perform a binary search to identify the interval containing the point to evaluate.
int lowerx = 0;
int upperx = xsize-1;
while (upperx-lowerx > 1) {
int middle = (upperx+lowerx)/2;
if (x[middle] > u)
upperx = middle;
else
lowerx = middle;
}
int lowery = 0;
int uppery = ysize-1;
while (uppery-lowery > 1) {
int middle = (uppery+lowery)/2;
if (y[middle] > v)
uppery = middle;
else
lowery = middle;
}
double deltax = x[upperx]-x[lowerx];
double deltay = y[uppery]-y[lowery];
double da = (u-x[lowerx])/deltax;
double db = (v-y[lowery])/deltay;
const vector<double>& coeff = c[lowerx+(xsize-1)*lowery];
// Evaluate the spline to determine the derivatives.
dx = 0;
dy = 0;
for (int i = 3; i >= 0; i--) {
dx = db*dx + (3.0*coeff[i+3*4]*da + 2.0*coeff[i+2*4])*da + coeff[i+1*4];
dy = da*dy + (3.0*coeff[i*4+3]*db + 2.0*coeff[i*4+2])*db + coeff[i*4+1];
}
dx /= deltax;
dy /= deltay;
}
void SplineFitter::create3DNaturalSpline(const vector<double>& x, const vector<double>& y, const vector<double>& z, const vector<double>& values, vector<vector<double> >& c) {
int xsize = x.size(), ysize = y.size(), zsize = z.size();
int xysize = xsize*ysize;
if (xsize < 2 || ysize < 2 || zsize < 2)
throw OpenMMException("create2DNaturalSpline: must have at least two points along each axis");
if (values.size() != xsize*ysize*zsize)
throw OpenMMException("create2DNaturalSpline: incorrect number of values");
vector<double> d1(xsize*ysize*zsize), d2(xsize*ysize*zsize), d3(xsize*ysize*zsize);
vector<double> d12(xsize*ysize*zsize), d13(xsize*ysize*zsize), d23(xsize*ysize*zsize), d123(xsize*ysize*zsize);
vector<double> t(xsize), deriv(xsize);
// Compute derivatives with respect to x.
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < zsize; j++) {
for (int k = 0; k < xsize; k++)
t[k] = values[k+xsize*i+xysize*j];
SplineFitter::createNaturalSpline(x, t, deriv);
for (int k = 0; k < xsize; k++)
d1[k+xsize*i+xysize*j] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[k]);
}
}
// Compute derivatives with respect to y.
t.resize(ysize);
deriv.resize(ysize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < zsize; j++) {
for (int k = 0; k < ysize; k++)
t[k] = values[i+xsize*k+xysize*j];
SplineFitter::createNaturalSpline(y, t, deriv);
for (int k = 0; k < ysize; k++)
d2[i+xsize*k+xysize*j] = SplineFitter::evaluateSplineDerivative(y, t, deriv, y[k]);
}
}
// Compute derivatives with respect to z.
t.resize(zsize);
deriv.resize(zsize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
for (int k = 0; k < zsize; k++)
t[k] = values[i+xsize*j+xysize*k];
SplineFitter::createNaturalSpline(z, t, deriv);
for (int k = 0; k < zsize; k++)
d3[i+xsize*j+xysize*k] = SplineFitter::evaluateSplineDerivative(z, t, deriv, z[k]);
}
}
// Compute second derivatives with respect to x and y.
t.resize(xsize);
deriv.resize(xsize);
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < zsize; j++) {
for (int k = 0; k < xsize; k++)
t[k] = d2[k+xsize*i+xysize*j];
SplineFitter::createNaturalSpline(x, t, deriv);
for (int k = 0; k < xsize; k++)
d12[k+xsize*i+xysize*j] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[k]);
}
}
// Compute second derivatives with respect to y and z.
t.resize(ysize);
deriv.resize(ysize);
for (int i = 0; i < zsize; i++) {
for (int j = 0; j < xsize; j++) {
for (int k = 0; k < ysize; k++)
t[k] = d3[j+xsize*k+xysize*i];
SplineFitter::createNaturalSpline(y, t, deriv);
for (int k = 0; k < ysize; k++)
d23[j+xsize*k+xysize*i] = SplineFitter::evaluateSplineDerivative(y, t, deriv, y[k]);
}
}
// Compute second derivatives with respect to x and z.
t.resize(zsize);
deriv.resize(zsize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++) {
for (int k = 0; k < zsize; k++)
t[k] = d1[i+xsize*j+xysize*k];
SplineFitter::createNaturalSpline(z, t, deriv);
for (int k = 0; k < zsize; k++)
d13[i+xsize*j+xysize*k] = SplineFitter::evaluateSplineDerivative(z, t, deriv, z[k]);
}
}
// Compute third derivatives with respect to x, y, and z.
t.resize(xsize);
deriv.resize(xsize);
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < zsize; j++) {
for (int k = 0; k < xsize; k++)
t[k] = d23[k+xsize*i+xysize*j];
SplineFitter::createNaturalSpline(x, t, deriv);
for (int k = 0; k < xsize; k++)
d123[k+xsize*i+xysize*j] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[k]);
}
}
// Now compute the coefficients. This involves multiplying by a sparse 64x64 matrix, given
// here in packed form.
const int wt[] = {
1,0,1,
1,8,1,
4,0,-3,1,3,8,-2,9,-1,
4,0,2,1,-2,8,1,9,1,
1,16,1,
1,32,1,
4,16,-3,17,3,32,-2,33,-1,
4,16,2,17,-2,32,1,33,1,
4,0,-3,2,3,16,-2,18,-1,
4,8,-3,10,3,32,-2,34,-1,
16,0,9,1,-9,2,-9,3,9,8,6,9,3,10,-6,11,-3,16,6,17,-6,18,3,19,-3,32,4,33,2,34,2,35,1,
16,0,-6,1,6,2,6,3,-6,8,-3,9,-3,10,3,11,3,16,-4,17,4,18,-2,19,2,32,-2,33,-2,34,-1,35,-1,
4,0,2,2,-2,16,1,18,1,
4,8,2,10,-2,32,1,34,1,
16,0,-6,1,6,2,6,3,-6,8,-4,9,-2,10,4,11,2,16,-3,17,3,18,-3,19,3,32,-2,33,-1,34,-2,35,-1,
16,0,4,1,-4,2,-4,3,4,8,2,9,2,10,-2,11,-2,16,2,17,-2,18,2,19,-2,32,1,33,1,34,1,35,1,
1,24,1,
1,40,1,
4,24,-3,25,3,40,-2,41,-1,
4,24,2,25,-2,40,1,41,1,
1,48,1,
1,56,1,
4,48,-3,49,3,56,-2,57,-1,
4,48,2,49,-2,56,1,57,1,
4,24,-3,26,3,48,-2,50,-1,
4,40,-3,42,3,56,-2,58,-1,
16,24,9,25,-9,26,-9,27,9,40,6,41,3,42,-6,43,-3,48,6,49,-6,50,3,51,-3,56,4,57,2,58,2,59,1,
16,24,-6,25,6,26,6,27,-6,40,-3,41,-3,42,3,43,3,48,-4,49,4,50,-2,51,2,56,-2,57,-2,58,-1,59,-1,
4,24,2,26,-2,48,1,50,1,
4,40,2,42,-2,56,1,58,1,
16,24,-6,25,6,26,6,27,-6,40,-4,41,-2,42,4,43,2,48,-3,49,3,50,-3,51,3,56,-2,57,-1,58,-2,59,-1,
16,24,4,25,-4,26,-4,27,4,40,2,41,2,42,-2,43,-2,48,2,49,-2,50,2,51,-2,56,1,57,1,58,1,59,1,
4,0,-3,4,3,24,-2,28,-1,
4,8,-3,12,3,40,-2,44,-1,
16,0,9,1,-9,4,-9,5,9,8,6,9,3,12,-6,13,-3,24,6,25,-6,28,3,29,-3,40,4,41,2,44,2,45,1,
16,0,-6,1,6,4,6,5,-6,8,-3,9,-3,12,3,13,3,24,-4,25,4,28,-2,29,2,40,-2,41,-2,44,-1,45,-1,
4,16,-3,20,3,48,-2,52,-1,
4,32,-3,36,3,56,-2,60,-1,
16,16,9,17,-9,20,-9,21,9,32,6,33,3,36,-6,37,-3,48,6,49,-6,52,3,53,-3,56,4,57,2,60,2,61,1,
16,16,-6,17,6,20,6,21,-6,32,-3,33,-3,36,3,37,3,48,-4,49,4,52,-2,53,2,56,-2,57,-2,60,-1,61,-1,
16,0,9,2,-9,4,-9,6,9,16,6,18,3,20,-6,22,-3,24,6,26,-6,28,3,30,-3,48,4,50,2,52,2,54,1,
16,8,9,10,-9,12,-9,14,9,32,6,34,3,36,-6,38,-3,40,6,42,-6,44,3,46,-3,56,4,58,2,60,2,62,1,
64,0,-27,1,27,2,27,3,-27,4,27,5,-27,6,-27,7,27,8,-18,9,-9,10,18,11,9,12,18,13,9,14,-18,15,-9,16,-18,17,18,18,-9,19,9,20,18,21,-18,22,9,23,-9,24,-18,25,18,26,18,27,-18,28,-9,29,9,30,9,31,-9,32,-12,33,-6,34,-6,35,-3,36,12,37,6,38,6,39,3,40,-12,41,-6,42,12,43,6,44,-6,45,-3,46,6,47,3,48,-12,49,12,50,-6,51,6,52,-6,53,6,54,-3,55,3,56,-8,57,-4,58,-4,59,-2,60,-4,61,-2,62,-2,63,-1,
64,0,18,1,-18,2,-18,3,18,4,-18,5,18,6,18,7,-18,8,9,9,9,10,-9,11,-9,12,-9,13,-9,14,9,15,9,16,12,17,-12,18,6,19,-6,20,-12,21,12,22,-6,23,6,24,12,25,-12,26,-12,27,12,28,6,29,-6,30,-6,31,6,32,6,33,6,34,3,35,3,36,-6,37,-6,38,-3,39,-3,40,6,41,6,42,-6,43,-6,44,3,45,3,46,-3,47,-3,48,8,49,-8,50,4,51,-4,52,4,53,-4,54,2,55,-2,56,4,57,4,58,2,59,2,60,2,61,2,62,1,63,1,
16,0,-6,2,6,4,6,6,-6,16,-3,18,-3,20,3,22,3,24,-4,26,4,28,-2,30,2,48,-2,50,-2,52,-1,54,-1,
16,8,-6,10,6,12,6,14,-6,32,-3,34,-3,36,3,38,3,40,-4,42,4,44,-2,46,2,56,-2,58,-2,60,-1,62,-1,
64,0,18,1,-18,2,-18,3,18,4,-18,5,18,6,18,7,-18,8,12,9,6,10,-12,11,-6,12,-12,13,-6,14,12,15,6,16,9,17,-9,18,9,19,-9,20,-9,21,9,22,-9,23,9,24,12,25,-12,26,-12,27,12,28,6,29,-6,30,-6,31,6,32,6,33,3,34,6,35,3,36,-6,37,-3,38,-6,39,-3,40,8,41,4,42,-8,43,-4,44,4,45,2,46,-4,47,-2,48,6,49,-6,50,6,51,-6,52,3,53,-3,54,3,55,-3,56,4,57,2,58,4,59,2,60,2,61,1,62,2,63,1,
64,0,-12,1,12,2,12,3,-12,4,12,5,-12,6,-12,7,12,8,-6,9,-6,10,6,11,6,12,6,13,6,14,-6,15,-6,16,-6,17,6,18,-6,19,6,20,6,21,-6,22,6,23,-6,24,-8,25,8,26,8,27,-8,28,-4,29,4,30,4,31,-4,32,-3,33,-3,34,-3,35,-3,36,3,37,3,38,3,39,3,40,-4,41,-4,42,4,43,4,44,-2,45,-2,46,2,47,2,48,-4,49,4,50,-4,51,4,52,-2,53,2,54,-2,55,2,56,-2,57,-2,58,-2,59,-2,60,-1,61,-1,62,-1,63,-1,
4,0,2,4,-2,24,1,28,1,
4,8,2,12,-2,40,1,44,1,
16,0,-6,1,6,4,6,5,-6,8,-4,9,-2,12,4,13,2,24,-3,25,3,28,-3,29,3,40,-2,41,-1,44,-2,45,-1,
16,0,4,1,-4,4,-4,5,4,8,2,9,2,12,-2,13,-2,24,2,25,-2,28,2,29,-2,40,1,41,1,44,1,45,1,
4,16,2,20,-2,48,1,52,1,
4,32,2,36,-2,56,1,60,1,
16,16,-6,17,6,20,6,21,-6,32,-4,33,-2,36,4,37,2,48,-3,49,3,52,-3,53,3,56,-2,57,-1,60,-2,61,-1,
16,16,4,17,-4,20,-4,21,4,32,2,33,2,36,-2,37,-2,48,2,49,-2,52,2,53,-2,56,1,57,1,60,1,61,1,
16,0,-6,2,6,4,6,6,-6,16,-4,18,-2,20,4,22,2,24,-3,26,3,28,-3,30,3,48,-2,50,-1,52,-2,54,-1,
16,8,-6,10,6,12,6,14,-6,32,-4,34,-2,36,4,38,2,40,-3,42,3,44,-3,46,3,56,-2,58,-1,60,-2,62,-1,
64,0,18,1,-18,2,-18,3,18,4,-18,5,18,6,18,7,-18,8,12,9,6,10,-12,11,-6,12,-12,13,-6,14,12,15,6,16,12,17,-12,18,6,19,-6,20,-12,21,12,22,-6,23,6,24,9,25,-9,26,-9,27,9,28,9,29,-9,30,-9,31,9,32,8,33,4,34,4,35,2,36,-8,37,-4,38,-4,39,-2,40,6,41,3,42,-6,43,-3,44,6,45,3,46,-6,47,-3,48,6,49,-6,50,3,51,-3,52,6,53,-6,54,3,55,-3,56,4,57,2,58,2,59,1,60,4,61,2,62,2,63,1,
64,0,-12,1,12,2,12,3,-12,4,12,5,-12,6,-12,7,12,8,-6,9,-6,10,6,11,6,12,6,13,6,14,-6,15,-6,16,-8,17,8,18,-4,19,4,20,8,21,-8,22,4,23,-4,24,-6,25,6,26,6,27,-6,28,-6,29,6,30,6,31,-6,32,-4,33,-4,34,-2,35,-2,36,4,37,4,38,2,39,2,40,-3,41,-3,42,3,43,3,44,-3,45,-3,46,3,47,3,48,-4,49,4,50,-2,51,2,52,-4,53,4,54,-2,55,2,56,-2,57,-2,58,-1,59,-1,60,-2,61,-2,62,-1,63,-1,
16,0,4,2,-4,4,-4,6,4,16,2,18,2,20,-2,22,-2,24,2,26,-2,28,2,30,-2,48,1,50,1,52,1,54,1,
16,8,4,10,-4,12,-4,14,4,32,2,34,2,36,-2,38,-2,40,2,42,-2,44,2,46,-2,56,1,58,1,60,1,62,1,
64,0,-12,1,12,2,12,3,-12,4,12,5,-12,6,-12,7,12,8,-8,9,-4,10,8,11,4,12,8,13,4,14,-8,15,-4,16,-6,17,6,18,-6,19,6,20,6,21,-6,22,6,23,-6,24,-6,25,6,26,6,27,-6,28,-6,29,6,30,6,31,-6,32,-4,33,-2,34,-4,35,-2,36,4,37,2,38,4,39,2,40,-4,41,-2,42,4,43,2,44,-4,45,-2,46,4,47,2,48,-3,49,3,50,-3,51,3,52,-3,53,3,54,-3,55,3,56,-2,57,-1,58,-2,59,-1,60,-2,61,-1,62,-2,63,-1,
64,0,8,1,-8,2,-8,3,8,4,-8,5,8,6,8,7,-8,8,4,9,4,10,-4,11,-4,12,-4,13,-4,14,4,15,4,16,4,17,-4,18,4,19,-4,20,-4,21,4,22,-4,23,4,24,4,25,-4,26,-4,27,4,28,4,29,-4,30,-4,31,4,32,2,33,2,34,2,35,2,36,-2,37,-2,38,-2,39,-2,40,2,41,2,42,-2,43,-2,44,2,45,2,46,-2,47,-2,48,2,49,-2,50,2,51,-2,52,2,53,-2,54,2,55,-2,56,1,57,1,58,1,59,1,60,1,61,1,62,1,63,1
};
vector<vector<int> > weight(64);
int index = 0;
for (int i = 0; i < 64; i++) {
int numElements = wt[index++];
for (int j = 0; j < numElements; j++) {
weight[i].push_back(wt[index++]);
weight[i].push_back(wt[index++]);
}
}
vector<double> rhs(64);
c.resize((xsize-1)*(ysize-1)*(zsize-1));
for (int i = 0; i < xsize-1; i++) {
for (int j = 0; j < ysize-1; j++) {
for (int k = 0; k < zsize-1; k++) {
// Compute the 64 coefficients for patch (i, j, k).
int nexti = i+1;
int nextj = j+1;
int nextk = k+1;
double deltax = x[nexti]-x[i];
double deltay = y[nextj]-y[j];
double deltaz = z[nextj]-z[j];
double e[] = {values[i+j*xsize+k*xysize], values[nexti+j*xsize+k*xysize], values[i+nextj*xsize+k*xysize], values[nexti+nextj*xsize+k*xysize], values[i+j*xsize+nextk*xysize], values[nexti+j*xsize+nextk*xysize], values[i+nextj*xsize+nextk*xysize], values[nexti+nextj*xsize+nextk*xysize]};
double e1[] = {d1[i+j*xsize+k*xysize], d1[nexti+j*xsize+k*xysize], d1[i+nextj*xsize+k*xysize], d1[nexti+nextj*xsize+k*xysize], d1[i+j*xsize+nextk*xysize], d1[nexti+j*xsize+nextk*xysize], d1[i+nextj*xsize+nextk*xysize], d1[nexti+nextj*xsize+nextk*xysize]};
double e2[] = {d2[i+j*xsize+k*xysize], d2[nexti+j*xsize+k*xysize], d2[i+nextj*xsize+k*xysize], d2[nexti+nextj*xsize+k*xysize], d2[i+j*xsize+nextk*xysize], d2[nexti+j*xsize+nextk*xysize], d2[i+nextj*xsize+nextk*xysize], d2[nexti+nextj*xsize+nextk*xysize]};
double e3[] = {d3[i+j*xsize+k*xysize], d3[nexti+j*xsize+k*xysize], d3[i+nextj*xsize+k*xysize], d3[nexti+nextj*xsize+k*xysize], d3[i+j*xsize+nextk*xysize], d3[nexti+j*xsize+nextk*xysize], d3[i+nextj*xsize+nextk*xysize], d3[nexti+nextj*xsize+nextk*xysize]};
double e12[] = {d12[i+j*xsize+k*xysize], d12[nexti+j*xsize+k*xysize], d12[i+nextj*xsize+k*xysize], d12[nexti+nextj*xsize+k*xysize], d12[i+j*xsize+nextk*xysize], d12[nexti+j*xsize+nextk*xysize], d12[i+nextj*xsize+nextk*xysize], d12[nexti+nextj*xsize+nextk*xysize]};
double e13[] = {d13[i+j*xsize+k*xysize], d13[nexti+j*xsize+k*xysize], d13[i+nextj*xsize+k*xysize], d13[nexti+nextj*xsize+k*xysize], d13[i+j*xsize+nextk*xysize], d13[nexti+j*xsize+nextk*xysize], d13[i+nextj*xsize+nextk*xysize], d13[nexti+nextj*xsize+nextk*xysize]};
double e23[] = {d23[i+j*xsize+k*xysize], d23[nexti+j*xsize+k*xysize], d23[i+nextj*xsize+k*xysize], d23[nexti+nextj*xsize+k*xysize], d23[i+j*xsize+nextk*xysize], d23[nexti+j*xsize+nextk*xysize], d23[i+nextj*xsize+nextk*xysize], d23[nexti+nextj*xsize+nextk*xysize]};
double e123[] = {d123[i+j*xsize+k*xysize], d123[nexti+j*xsize+k*xysize], d123[i+nextj*xsize+k*xysize], d123[nexti+nextj*xsize+k*xysize], d123[i+j*xsize+nextk*xysize], d123[nexti+j*xsize+nextk*xysize], d123[i+nextj*xsize+nextk*xysize], d123[nexti+nextj*xsize+nextk*xysize]};
for (int m = 0; m < 8; m++) {
rhs[m] = e[m];
rhs[m+8] = e1[m]*deltax;
rhs[m+16] = e2[m]*deltay;
rhs[m+24] = e3[m]*deltaz;
rhs[m+32] = e12[m]*deltax*deltay;
rhs[m+40] = e13[m]*deltax*deltaz;
rhs[m+48] = e23[m]*deltay*deltaz;
rhs[m+56] = e123[m]*deltax*deltay*deltaz;
}
vector<double>& coeff = c[i+j*(xsize-1)+k*(xsize-1)*(ysize-1)];
coeff.resize(64);
for (int m = 0; m < 64; m++) {
double sum = 0.0;
int numElements = weight[m].size();
for (int n = 0; n < numElements; n += 2)
sum += weight[m][n+1]*rhs[weight[m][n]];
coeff[m] = sum;
}
}
}
}
}
double SplineFitter::evaluate3DSpline(const vector<double>& x, const vector<double>& y, const vector<double>& z, const vector<double>& values, const vector<vector<double> >& c, double u, double v, double w) {
int xsize = x.size();
int ysize = y.size();
int zsize = z.size();
if (u < x[0] || u > x[xsize-1] || v < y[0] || v > y[ysize-1] || w < z[0] || w > z[zsize-1])
throw OpenMMException("evaluate3DSpline: specified point is outside the range defined by the spline");
// Perform a binary search to identify the interval containing the point to evaluate.
int lowerx = 0;
int upperx = xsize-1;
while (upperx-lowerx > 1) {
int middle = (upperx+lowerx)/2;
if (x[middle] > u)
upperx = middle;
else
lowerx = middle;
}
int lowery = 0;
int uppery = ysize-1;
while (uppery-lowery > 1) {
int middle = (uppery+lowery)/2;
if (y[middle] > v)
uppery = middle;
else
lowery = middle;
}
int lowerz = 0;
int upperz = zsize-1;
while (upperz-lowerz > 1) {
int middle = (upperz+lowerz)/2;
if (z[middle] > w)
upperz = middle;
else
lowerz = middle;
}
double deltax = x[upperx]-x[lowerx];
double deltay = y[uppery]-y[lowery];
double deltaz = z[upperz]-z[lowerz];
double da = (u-x[lowerx])/deltax;
double db = (v-y[lowery])/deltay;
double dc = (w-z[lowerz])/deltaz;
const vector<double>& coeff = c[lowerx+(xsize-1)*lowery+(xsize-1)*(ysize-1)*lowerz];
// Evaluate the spline to determine the value and gradients.
double value[] = {0, 0, 0, 0};
for (int i = 3; i >= 0; i--) {
for (int j = 0; j < 4; j++) {
int base = 4*i + 16*j;
value[j] = db*value[j] + ((coeff[base+3]*da + coeff[base+2])*da + coeff[base+1])*da + coeff[base];
}
}
return value[0] + dc*(value[1] + dc*(value[2] + dc*value[3]));
}
void SplineFitter::evaluate3DSplineDerivatives(const vector<double>& x, const vector<double>& y, const vector<double>& z, const vector<double>& values, const vector<vector<double> >& c, double u, double v, double w, double& dx, double& dy, double& dz) {
int xsize = x.size();
int ysize = y.size();
int zsize = z.size();
if (u < x[0] || u > x[xsize-1] || v < y[0] || v > y[ysize-1] || w < z[0] || w > z[zsize-1])
throw OpenMMException("evaluate3DSpline: specified point is outside the range defined by the spline");
// Perform a binary search to identify the interval containing the point to evaluate.
int lowerx = 0;
int upperx = xsize-1;
while (upperx-lowerx > 1) {
int middle = (upperx+lowerx)/2;
if (x[middle] > u)
upperx = middle;
else
lowerx = middle;
}
int lowery = 0;
int uppery = ysize-1;
while (uppery-lowery > 1) {
int middle = (uppery+lowery)/2;
if (y[middle] > v)
uppery = middle;
else
lowery = middle;
}
int lowerz = 0;
int upperz = zsize-1;
while (upperz-lowerz > 1) {
int middle = (upperz+lowerz)/2;
if (z[middle] > w)
upperz = middle;
else
lowerz = middle;
}
double deltax = x[upperx]-x[lowerx];
double deltay = y[uppery]-y[lowery];
double deltaz = z[upperz]-z[lowerz];
double da = (u-x[lowerx])/deltax;
double db = (v-y[lowery])/deltay;
double dc = (w-z[lowerz])/deltaz;
const vector<double>& coeff = c[lowerx+(xsize-1)*lowery+(xsize-1)*(ysize-1)*lowerz];
// Evaluate the spline to determine the derivatives.
double derivx[] = {0, 0, 0, 0};
double derivy[] = {0, 0, 0, 0};
double derivz[] = {0, 0, 0, 0};
for (int i = 3; i >= 0; i--) {
for (int j = 0; j < 4; j++) {
int base = 4*i + 16*j;
derivx[j] = db*derivx[j] + (3.0*coeff[base+3]*da + 2.0*coeff[base+2])*da + coeff[base+1];
derivz[j] = db*derivz[j] + ((coeff[base+3]*da + coeff[base+2])*da + coeff[base+1])*da + coeff[base];
base = i + 16*j;
derivy[j] = da*derivy[j] + (3.0*coeff[base+12]*db + 2.0*coeff[base+8])*db + coeff[base+4];
}
}
dx = derivx[0] + dc*(derivx[1] + dc*(derivx[2] + dc*derivx[3]));
dy = derivy[0] + dc*(derivy[1] + dc*(derivy[2] + dc*derivy[3]));
dz = derivz[1] + dc*(2.0*derivz[2] + 3.0*dc*derivz[3]);
dx /= deltax;
dy /= deltay;
dz /= deltaz;
}
openmmapi/src/TabulatedFunction.cpp 0 → 100644
View file @ eb9f735a
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2014 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/TabulatedFunction.h"
#include "openmm/OpenMMException.h"
using namespace OpenMM;
using namespace std;
Continuous1DFunction::Continuous1DFunction(const vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("Continuous1DFunction: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("Continuous1DFunction: a tabulated function must have at least two points");
this->values = values;
this->min = min;
this->max = max;
}
void Continuous1DFunction::getFunctionParameters(vector<double>& values, double& min, double& max) const {
values = this->values;
min = this->min;
max = this->max;
}
void Continuous1DFunction::setFunctionParameters(const vector<double>& values, double min, double max) {
if (max <= min)
throw OpenMMException("Continuous1DFunction: max <= min for a tabulated function.");
if (values.size() < 2)
throw OpenMMException("Continuous1DFunction: a tabulated function must have at least two points");
this->values = values;
this->min = min;
this->max = max;
}
Continuous2DFunction::Continuous2DFunction(int xsize, int ysize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax) {
if (xsize < 2 || ysize < 2)
throw OpenMMException("Continuous2DFunction: must have at least two points along each axis");
if (values.size() != xsize*ysize)
throw OpenMMException("Continuous2DFunction: incorrect number of values");
if (xmax <= xmin)
throw OpenMMException("Continuous2DFunction: xmax <= xmin for a tabulated function.");
if (ymax <= ymin)
throw OpenMMException("Continuous2DFunction: ymax <= ymin for a tabulated function.");
this->values = values;
this->xsize = xsize;
this->ysize = ysize;
this->xmin = xmin;
this->xmax = xmax;
this->ymin = ymin;
this->ymax = ymax;
}
void Continuous2DFunction::getFunctionParameters(int& xsize, int& ysize, vector<double>& values, double& xmin, double& xmax, double& ymin, double& ymax) const {
values = this->values;
xsize = this->xsize;
ysize = this->ysize;
xmin = this->xmin;
xmax = this->xmax;
ymin = this->ymin;
ymax = this->ymax;
}
void Continuous2DFunction::setFunctionParameters(int xsize, int ysize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax) {
if (xsize < 2 || ysize < 2)
throw OpenMMException("Continuous2DFunction: must have at least two points along each axis");
if (values.size() != xsize*ysize)
throw OpenMMException("Continuous2DFunction: incorrect number of values");
if (xmax <= xmin)
throw OpenMMException("Continuous2DFunction: xmax <= xmin for a tabulated function.");
if (ymax <= ymin)
throw OpenMMException("Continuous2DFunction: ymax <= ymin for a tabulated function.");
this->values = values;
this->xsize = xsize;
this->ysize = ysize;
this->xmin = xmin;
this->xmax = xmax;
this->ymin = ymin;
this->ymax = ymax;
}
Continuous3DFunction::Continuous3DFunction(int xsize, int ysize, int zsize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax) {
if (xsize < 2 || ysize < 2 || zsize < 2)
throw OpenMMException("Continuous3DFunction: must have at least two points along each axis");
if (values.size() != xsize*ysize*zsize)
throw OpenMMException("Continuous3DFunction: incorrect number of values");
if (xmax <= xmin)
throw OpenMMException("Continuous3DFunction: xmax <= xmin for a tabulated function.");
if (ymax <= ymin)
throw OpenMMException("Continuous3DFunction: ymax <= ymin for a tabulated function.");
if (zmax <= zmin)
throw OpenMMException("Continuous3DFunction: zmax <= zmin for a tabulated function.");
this->values = values;
this->xsize = xsize;
this->ysize = ysize;
this->zsize = zsize;
this->xmin = xmin;
this->xmax = xmax;
this->ymin = ymin;
this->ymax = ymax;
this->zmin = zmin;
this->zmax = zmax;
}
void Continuous3DFunction::getFunctionParameters(int& xsize, int& ysize, int& zsize, vector<double>& values, double& xmin, double& xmax, double& ymin, double& ymax, double& zmin, double& zmax) const {
values = this->values;
xsize = this->xsize;
ysize = this->ysize;
zsize = this->zsize;
xmin = this->xmin;
xmax = this->xmax;
ymin = this->ymin;
ymax = this->ymax;
zmin = this->zmin;
zmax = this->zmax;
}
void Continuous3DFunction::setFunctionParameters(int xsize, int ysize, int zsize, const vector<double>& values, double xmin, double xmax, double ymin, double ymax, double zmin, double zmax) {
if (xsize < 2 || ysize < 2 || zsize < 2)
throw OpenMMException("Continuous3DFunction: must have at least two points along each axis");
if (values.size() != xsize*ysize*zsize)
throw OpenMMException("Continuous3DFunction: incorrect number of values");
if (xmax <= xmin)
throw OpenMMException("Continuous3DFunction: xmax <= xmin for a tabulated function.");
if (ymax <= ymin)
throw OpenMMException("Continuous3DFunction: ymax <= ymin for a tabulated function.");
if (zmax <= zmin)
throw OpenMMException("Continuous3DFunction: zmax <= zmin for a tabulated function.");
this->values = values;
this->xsize = xsize;
this->ysize = ysize;
this->zsize = zsize;
this->xmin = xmin;
this->xmax = xmax;
this->ymin = ymin;
this->ymax = ymax;
this->zmin = zmin;
this->zmax = zmax;
}
Discrete1DFunction::Discrete1DFunction(const vector<double>& values) {
this->values = values;
}
void Discrete1DFunction::getFunctionParameters(vector<double>& values) const {
values = this->values;
}
void Discrete1DFunction::setFunctionParameters(const vector<double>& values) {
this->values = values;
}
Discrete2DFunction::Discrete2DFunction(int xsize, int ysize, const vector<double>& values) {
if (values.size() != xsize*ysize)
throw OpenMMException("Discrete2DFunction: incorrect number of values");
this->xsize = xsize;
this->ysize = ysize;
this->values = values;
}
void Discrete2DFunction::getFunctionParameters(int& xsize, int& ysize, vector<double>& values) const {
xsize = this->xsize;
ysize = this->ysize;
values = this->values;
}
void Discrete2DFunction::setFunctionParameters(int xsize, int ysize, const vector<double>& values) {
if (values.size() != xsize*ysize)
throw OpenMMException("Discrete2DFunction: incorrect number of values");
this->xsize = xsize;
this->ysize = ysize;
this->values = values;
}
Discrete3DFunction::Discrete3DFunction(int xsize, int ysize, int zsize, const vector<double>& values) {
if (values.size() != xsize*ysize*zsize)
throw OpenMMException("Discrete3DFunction: incorrect number of values");
this->xsize = xsize;
this->ysize = ysize;
this->zsize = zsize;
this->values = values;
}
void Discrete3DFunction::getFunctionParameters(int& xsize, int& ysize, int& zsize, vector<double>& values) const {
xsize = this->xsize;
ysize = this->ysize;
zsize = this->zsize;
values = this->values;
}
void Discrete3DFunction::setFunctionParameters(int xsize, int ysize, int zsize, const vector<double>& values) {
if (values.size() != xsize*ysize*zsize)
throw OpenMMException("Discrete3DFunction: incorrect number of values");
this->xsize = xsize;
this->ysize = ysize;
this->zsize = zsize;
this->values = values;
}
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