Skip to content

GitLab

"docs-source/vscode:/vscode.git/clone" did not exist on "de112a824e9064845637bcd05680e6cca828efed"
Commit 0827ad50 authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #1858 from peastman/localcoordinates 

LocalCoordinatesSite can depend on arbitary number of particles
parents 15b9fb48 ea7f44b2
Show whitespace changes
Inline Side-by-side
Showing with 564 additions and 419 deletions
docs-source/usersguide/application.rst
View file @ 0827ad50
......@@ -2051,10 +2051,9 @@ be represented with a virtual site. The :code:`type` attribute may equal
:code:`"localCoords"`\ , which correspond to the TwoParticleAverageSite, ThreeParticleAverageSite,
OutOfPlaneSite, and LocalCoordinatesSite classes respectively. The :code:`siteName`
attribute gives the name of the atom to represent with a virtual site. The atoms
it is calculated based on are specified by :code:`atomName1`\ , :code:`atomName2`\ ,
and (for virtual site classes that involve three atoms) :code:`atomName3`\ .
it is calculated based on are specified by :code:`atomName1`\ , :code:`atomName2`\ , etc.
(Some old force fields use the deprecated tags :code:`index`, :code:`atom1`,
:code:`atom2`, and :code:`atom3` to refer to them by index instead of name.)
:code:`atom2`, etc. to refer to them by index instead of name.)
The remaining attributes are specific to the virtual site class, and specify the
parameters for calculating the site position. For a TwoParticleAverageSite,
......@@ -2062,9 +2061,10 @@ they are :code:`weight1` and :code:`weight2`\ . For a
ThreeParticleAverageSite, they are :code:`weight1`\ , :code:`weight2`\ , and
\ :code:`weight3`\ . For an OutOfPlaneSite, they are :code:`weight12`\ ,
:code:`weight13`\ , and :code:`weightCross`\ . For a LocalCoordinatesSite, they
are :code:`wo1`\ , :code:`wo2`\ , :code:`wo3`\ , :code:`wx1`\ , :code:`wx2`\ ,
:code:`wx3`\ , :code:`wy1`\ , :code:`wy2`\ , :code:`wy3`\ , :code:`p1`\ ,
:code:`p2`\ , and :code:`p3`\ .
are :code:`p1`\ , :code:`p2`\ , and :code:`p3` (giving the x, y, and z coordinates
of the site position in the local coordinate system), and :code:`wo1`\ ,
:code:`wx1`\ , :code:`wy1`\ , :code:`wo2`\ , :code:`wx2`\ , :code:`wy2`\ , ...
(giving the weights for computing the origin, x axis, and y axis).
<Patches>
=========
......
docs-source/usersguide/theory.rst
View file @ 0827ad50
......@@ -1572,17 +1572,18 @@ specific types of rules. They are:
:math:`\mathbf{r}_{13} = \mathbf{r}_{3}-\mathbf{r}_{1}`\ . This allows
the virtual site to be located outside the plane of the three particles.
* LocalCoordinatesSite: The locations of three other particles are used to compute a local
* LocalCoordinatesSite: The locations of several other particles are used to compute a local
coordinate system, and the virtual site is placed at a fixed location in that coordinate
system. The origin of the coordinate system and the directions of its x and y axes
are each specified as a weighted sum of the locations of the three particles:
system. The number of particles used to define the coordinate system is user defined.
The origin of the coordinate system and the directions of its x and y axes
are each specified as a weighted sum of the locations of the other particles:
.. math::
\mathbf{o}={w}^{o}_{1}\mathbf{r}_{1} + {w}^{o}_{2}\mathbf{r}_{2} + {w}^{o}_{3}\mathbf{r}_{3}
\mathbf{o}={w}^{o}_{1}\mathbf{r}_{1} + {w}^{o}_{2}\mathbf{r}_{2} + ...
\mathbf{dx}={w}^{x}_{1}\mathbf{r}_{1} + {w}^{x}_{2}\mathbf{r}_{2} + {w}^{x}_{3}\mathbf{r}_{3}
\mathbf{dx}={w}^{x}_{1}\mathbf{r}_{1} + {w}^{x}_{2}\mathbf{r}_{2} + ...
\mathbf{dy}={w}^{y}_{1}\mathbf{r}_{1} + {w}^{y}_{2}\mathbf{r}_{2} + {w}^{y}_{3}\mathbf{r}_{3}
\mathbf{dy}={w}^{y}_{1}\mathbf{r}_{1} + {w}^{y}_{2}\mathbf{r}_{2} + ...
\mathbf{dz}=\mathbf{dx}\times \mathbf{dy}
..
......
openmmapi/include/openmm/VirtualSite.h
View file @ 0827ad50
......@@ -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) 2012-2014 Stanford University and the Authors. *
* Portions copyright (c) 2012-2017 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -168,21 +168,21 @@ private:
};
/**
* This is a VirtualSite that uses the locations of three other particles to compute a local
* This is a VirtualSite that uses the locations of several other particles to compute a local
* coordinate system, then places the virtual site at a fixed location in that coordinate
* system. The origin of the coordinate system and the directions of its x and y axes
* are each specified as a weighted sum of the locations of the three particles:
* are each specified as a weighted sum of the locations of the other particles:
*
* origin = w<sup>o</sup><sub>1</sub>r<sub>1</sub> + w<sup>o</sup><sub>2</sub>r<sub>2</sub> + w<sup>o</sup><sub>3</sub>r<sub>3</sub>
* origin = w<sup>o</sup><sub>1</sub>r<sub>1</sub> + w<sup>o</sup><sub>2</sub>r<sub>2</sub> + w<sup>o</sup><sub>3</sub>r<sub>3</sub> + ...
*
* xdir = w<sup>x</sup><sub>1</sub>r<sub>1</sub> + w<sup>x</sup><sub>2</sub>r<sub>2</sub> + w<sup>x</sup><sub>3</sub>r<sub>3</sub>
* xdir = w<sup>x</sup><sub>1</sub>r<sub>1</sub> + w<sup>x</sup><sub>2</sub>r<sub>2</sub> + w<sup>x</sup><sub>3</sub>r<sub>3</sub> + ...
*
* ydir = w<sup>y</sup><sub>1</sub>r<sub>1</sub> + w<sup>y</sup><sub>2</sub>r<sub>2</sub> + w<sup>y</sup><sub>3</sub>r<sub>3</sub>
* ydir = w<sup>y</sup><sub>1</sub>r<sub>1</sub> + w<sup>y</sup><sub>2</sub>r<sub>2</sub> + w<sup>y</sup><sub>3</sub>r<sub>3</sub> + ...
*
* For the origin, the three weights must add to one. For example if
* For the origin, the weights must add to one. For example if
* (w<sup>o</sup><sub>1</sub>, w<sup>o</sup><sub>2</sub>, w<sup>o</sup><sub>3</sub>) = (1.0, 0.0, 0.0),
* the origin of the local coordinate system is at the location of particle 1. For xdir and ydir,
* the weights must add to zero. For excample, if
* the weights must add to zero. For example, if
* (w<sup>x</sup><sub>1</sub>, w<sup>x</sup><sub>2</sub>, w<sup>x</sup><sub>3</sub>) = (-1.0, 0.5, 0.5),
* the x axis points from particle 1 toward the midpoint between particles 2 and 3.
*
......@@ -197,6 +197,17 @@ public:
/**
* Create a new LocalCoordinatesSite virtual site.
*
* @param particles the indices of the particle the site depends on
* @param originWeights the weight factors for the particles when computing the origin location
* @param xWeights the weight factors for the particles when computing xdir
* @param yWeights the weight factors for the particles when computing ydir
* @param localPosition the position of the virtual site in the local coordinate system
*/
LocalCoordinatesSite(const std::vector<int>& particles, const std::vector<double>& originWeights, const std::vector<double>& xWeights, const std::vector<double>& yWeights, const Vec3& localPosition);
/**
* Create a new LocalCoordinatesSite virtual site. This constructor assumes the site depends on
* exactly three other particles.
*
* @param particle1 the index of the first particle
* @param particle2 the index of the second particle
* @param particle3 the index of the third particle
......@@ -207,23 +218,42 @@ public:
*/
LocalCoordinatesSite(int particle1, int particle2, int particle3, const Vec3& originWeights, const Vec3& xWeights, const Vec3& yWeights, const Vec3& localPosition);
/**
* Get the weight factors for the three particles when computing the origin location.
* Get the weight factors for the particles when computing the origin location.
*/
void getOriginWeights(std::vector<double>& weights) const;
/**
* Get the weight factors for the particles when computing the origin location.
* This version assumes the site depends on exactly three other particles, and
* throws an exception if that is not the case.
*/
Vec3 getOriginWeights() const;
/**
* Get the weight factors for the particles when computing xdir.
*/
void getXWeights(std::vector<double>& weights) const;
/**
* Get the weight factors for the particles when computing xdir.
* This version assumes the site depends on exactly three other particles, and
* throws an exception if that is not the case.
*/
const Vec3& getOriginWeights() const;
Vec3 getXWeights() const;
/**
* Get the weight factors for the three particles when computing xdir.
* Get the weight factors for the particles when computing ydir.
*/
const Vec3& getXWeights() const;
void getYWeights(std::vector<double>& weights) const;
/**
* Get the weight factors for the three particles when computing ydir.
* Get the weight factors for the particles when computing ydir.
* This version assumes the site depends on exactly three other particles, and
* throws an exception if that is not the case.
*/
const Vec3& getYWeights() const;
Vec3 getYWeights() const;
/**
* Get the position of the virtual site in the local coordinate system.
*/
const Vec3& getLocalPosition() const;
private:
Vec3 originWeights, xWeights, yWeights, localPosition;
std::vector<double> originWeights, xWeights, yWeights;
Vec3 localPosition;
};
} // namespace OpenMM
......
openmmapi/src/VirtualSite.cpp
View file @ 0827ad50
......@@ -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) 2012-2014 Stanford University and the Authors. *
* Portions copyright (c) 2012-2017 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -105,8 +105,30 @@ double OutOfPlaneSite::getWeightCross() const {
return weightCross;
}
LocalCoordinatesSite::LocalCoordinatesSite(int particle1, int particle2, int particle3, const Vec3& originWeights, const Vec3& xWeights, const Vec3& yWeights, const Vec3& localPosition) :
LocalCoordinatesSite::LocalCoordinatesSite(const vector<int>& particles, const vector<double>& originWeights, const vector<double>& xWeights, const vector<double>& yWeights, const Vec3& localPosition) :
originWeights(originWeights), xWeights(xWeights), yWeights(yWeights), localPosition(localPosition) {
int numParticles = particles.size();
if (numParticles < 2)
throw OpenMMException("LocalCoordinatesSite: Must depend on at least two other particles");
if (originWeights.size() != numParticles || xWeights.size() != numParticles || yWeights.size() != numParticles)
throw OpenMMException("LocalCoordinatesSite: Number of weights does not match number of particles");
double originSum = 0, xSum = 0, ySum = 0;
for (int i = 0; i < numParticles; i++) {
originSum += originWeights[i];
xSum += xWeights[i];
ySum += yWeights[i];
}
if (fabs(originSum-1.0) > 1e-6)
throw OpenMMException("LocalCoordinatesSite: Weights for computing origin must add to 1");
if (fabs(xSum) > 1e-6)
throw OpenMMException("LocalCoordinatesSite: Weights for computing x axis must add to 0");
if (fabs(ySum) > 1e-6)
throw OpenMMException("LocalCoordinatesSite: Weights for computing y axis must add to 0");
setParticles(particles);
}
LocalCoordinatesSite::LocalCoordinatesSite(int particle1, int particle2, int particle3, const Vec3& originWeights, const Vec3& xWeights, const Vec3& yWeights, const Vec3& localPosition) :
localPosition(localPosition) {
if (fabs(originWeights[0]+originWeights[1]+originWeights[2]-1.0) > 1e-6)
throw OpenMMException("LocalCoordinatesSite: Weights for computing origin must add to 1");
if (fabs(xWeights[0]+xWeights[1]+xWeights[2]) > 1e-6)
......@@ -118,18 +140,45 @@ LocalCoordinatesSite::LocalCoordinatesSite(int particle1, int particle2, int par
particles[1] = particle2;
particles[2] = particle3;
setParticles(particles);
this->originWeights.push_back(originWeights[0]);
this->originWeights.push_back(originWeights[1]);
this->originWeights.push_back(originWeights[2]);
this->xWeights.push_back(xWeights[0]);
this->xWeights.push_back(xWeights[1]);
this->xWeights.push_back(xWeights[2]);
this->yWeights.push_back(yWeights[0]);
this->yWeights.push_back(yWeights[1]);
this->yWeights.push_back(yWeights[2]);
}
void LocalCoordinatesSite::getOriginWeights(vector<double>& weights) const {
weights = originWeights;
}
Vec3 LocalCoordinatesSite::getOriginWeights() const {
if (originWeights.size() != 3)
throw OpenMMException("LocalCoordinatesSite: This version of getOriginWeights() requires the site to depend on three particles");
return Vec3(originWeights[0], originWeights[1], originWeights[2]);
}
void LocalCoordinatesSite::getXWeights(vector<double>& weights) const {
weights = xWeights;
}
const Vec3& LocalCoordinatesSite::getOriginWeights() const {
return originWeights;
Vec3 LocalCoordinatesSite::getXWeights() const {
if (xWeights.size() != 3)
throw OpenMMException("LocalCoordinatesSite: This version of getXWeights() requires the site to depend on three particles");
return Vec3(xWeights[0], xWeights[1], xWeights[2]);
}
const Vec3& LocalCoordinatesSite::getXWeights() const {
return xWeights;
void LocalCoordinatesSite::getYWeights(vector<double>& weights) const {
weights = yWeights;
}
const Vec3& LocalCoordinatesSite::getYWeights() const {
return yWeights;
Vec3 LocalCoordinatesSite::getYWeights() const {
if (yWeights.size() != 3)
throw OpenMMException("LocalCoordinatesSite: This version of getYWeights() requires the site to depend on three particles");
return Vec3(yWeights[0], yWeights[1], yWeights[2]);
}
const Vec3& LocalCoordinatesSite::getLocalPosition() const {
......
platforms/cuda/include/CudaIntegrationUtilities.h
View file @ 0827ad50
......@@ -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) 2009-2014 Stanford University and the Authors. *
* Portions copyright (c) 2009-2017 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -158,8 +158,11 @@ private:
CudaArray* vsite3AvgWeights;
CudaArray* vsiteOutOfPlaneAtoms;
CudaArray* vsiteOutOfPlaneWeights;
CudaArray* vsiteLocalCoordsIndex;
CudaArray* vsiteLocalCoordsAtoms;
CudaArray* vsiteLocalCoordsParams;
CudaArray* vsiteLocalCoordsWeights;
CudaArray* vsiteLocalCoordsPos;
CudaArray* vsiteLocalCoordsStartIndex;
int randomPos;
int lastSeed, numVsites;
double2 lastStepSize;
......
platforms/cuda/src/CudaIntegrationUtilities.cpp
View file @ 0827ad50
......@@ -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) 2009-2015 Stanford University and the Authors. *
* Portions copyright (c) 2009-2017 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -103,7 +103,8 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
ccmaReducedMass(NULL), ccmaAtomConstraints(NULL), ccmaNumAtomConstraints(NULL), ccmaConstraintMatrixColumn(NULL),
ccmaConstraintMatrixValue(NULL), ccmaDelta1(NULL), ccmaDelta2(NULL), ccmaConverged(NULL), ccmaConvergedMemory(NULL),
vsite2AvgAtoms(NULL), vsite2AvgWeights(NULL), vsite3AvgAtoms(NULL), vsite3AvgWeights(NULL),
vsiteOutOfPlaneAtoms(NULL), vsiteOutOfPlaneWeights(NULL), vsiteLocalCoordsAtoms(NULL), vsiteLocalCoordsParams(NULL) {
vsiteOutOfPlaneAtoms(NULL), vsiteOutOfPlaneWeights(NULL), vsiteLocalCoordsIndex(NULL), vsiteLocalCoordsAtoms(NULL),
vsiteLocalCoordsWeights(NULL), vsiteLocalCoordsPos(NULL), vsiteLocalCoordsStartIndex(NULL) {
// Create workspace arrays.
lastStepSize = make_double2(0.0, 0.0);
......@@ -454,8 +455,11 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
vector<double4> vsite3AvgWeightVec;
vector<int4> vsiteOutOfPlaneAtomVec;
vector<double4> vsiteOutOfPlaneWeightVec;
vector<int4> vsiteLocalCoordsAtomVec;
vector<double> vsiteLocalCoordsParamVec;
vector<int> vsiteLocalCoordsIndexVec;
vector<int> vsiteLocalCoordsAtomVec;
vector<int> vsiteLocalCoordsStartVec;
vector<double> vsiteLocalCoordsWeightVec;
vector<double4> vsiteLocalCoordsPosVec;
for (int i = 0; i < numAtoms; i++) {
if (system.isVirtualSite(i)) {
if (dynamic_cast<const TwoParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
......@@ -480,64 +484,72 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
vsiteOutOfPlaneWeightVec.push_back(make_double4(site.getWeight12(), site.getWeight13(), site.getWeightCross(), 0.0));
}
else if (dynamic_cast<const LocalCoordinatesSite*>(&system.getVirtualSite(i)) != NULL) {
// An out of plane site.
// A local coordinates site.
const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
vsiteLocalCoordsAtomVec.push_back(make_int4(i, site.getParticle(0), site.getParticle(1), site.getParticle(2)));
Vec3 origin = site.getOriginWeights();
Vec3 x = site.getXWeights();
Vec3 y = site.getYWeights();
int numParticles = site.getNumParticles();
vector<double> origin, x, y;
site.getOriginWeights(origin);
site.getXWeights(x);
site.getYWeights(y);
vsiteLocalCoordsIndexVec.push_back(i);
vsiteLocalCoordsStartVec.push_back(vsiteLocalCoordsAtomVec.size());
for (int j = 0; j < numParticles; j++) {
vsiteLocalCoordsAtomVec.push_back(site.getParticle(j));
vsiteLocalCoordsWeightVec.push_back(origin[j]);
vsiteLocalCoordsWeightVec.push_back(x[j]);
vsiteLocalCoordsWeightVec.push_back(y[j]);
}
Vec3 pos = site.getLocalPosition();
vsiteLocalCoordsParamVec.push_back(origin[0]);
vsiteLocalCoordsParamVec.push_back(origin[1]);
vsiteLocalCoordsParamVec.push_back(origin[2]);
vsiteLocalCoordsParamVec.push_back(x[0]);
vsiteLocalCoordsParamVec.push_back(x[1]);
vsiteLocalCoordsParamVec.push_back(x[2]);
vsiteLocalCoordsParamVec.push_back(y[0]);
vsiteLocalCoordsParamVec.push_back(y[1]);
vsiteLocalCoordsParamVec.push_back(y[2]);
vsiteLocalCoordsParamVec.push_back(pos[0]);
vsiteLocalCoordsParamVec.push_back(pos[1]);
vsiteLocalCoordsParamVec.push_back(pos[2]);
vsiteLocalCoordsPosVec.push_back(make_double4(pos[0], pos[1], pos[2], 0.0));
}
}
}
vsiteLocalCoordsStartVec.push_back(vsiteLocalCoordsAtomVec.size());
int num2Avg = vsite2AvgAtomVec.size();
int num3Avg = vsite3AvgAtomVec.size();
int numOutOfPlane = vsiteOutOfPlaneAtomVec.size();
int numLocalCoords = vsiteLocalCoordsAtomVec.size();
int numLocalCoords = vsiteLocalCoordsPosVec.size();
vsite2AvgAtoms = CudaArray::create<int4>(context, max(1, num2Avg), "vsite2AvgAtoms");
vsite3AvgAtoms = CudaArray::create<int4>(context, max(1, num3Avg), "vsite3AvgAtoms");
vsiteOutOfPlaneAtoms = CudaArray::create<int4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneAtoms");
vsiteLocalCoordsAtoms = CudaArray::create<int4>(context, max(1, numLocalCoords), "vsiteLocalCoordinatesAtoms");
vsiteLocalCoordsIndex = CudaArray::create<int>(context, max(1, (int) vsiteLocalCoordsIndexVec.size()), "vsiteLocalCoordsIndex");
vsiteLocalCoordsAtoms = CudaArray::create<int>(context, max(1, (int) vsiteLocalCoordsAtomVec.size()), "vsiteLocalCoordsAtoms");
vsiteLocalCoordsStartIndex = CudaArray::create<int>(context, max(1, (int) vsiteLocalCoordsStartVec.size()), "vsiteLocalCoordsStartIndex");
if (num2Avg > 0)
vsite2AvgAtoms->upload(vsite2AvgAtomVec);
if (num3Avg > 0)
vsite3AvgAtoms->upload(vsite3AvgAtomVec);
if (numOutOfPlane > 0)
vsiteOutOfPlaneAtoms->upload(vsiteOutOfPlaneAtomVec);
if (numLocalCoords > 0)
if (numLocalCoords > 0) {
vsiteLocalCoordsIndex->upload(vsiteLocalCoordsIndexVec);
vsiteLocalCoordsAtoms->upload(vsiteLocalCoordsAtomVec);
vsiteLocalCoordsStartIndex->upload(vsiteLocalCoordsStartVec);
}
if (context.getUseDoublePrecision()) {
vsite2AvgWeights = CudaArray::create<double2>(context, max(1, num2Avg), "vsite2AvgWeights");
vsite3AvgWeights = CudaArray::create<double4>(context, max(1, num3Avg), "vsite3AvgWeights");
vsiteOutOfPlaneWeights = CudaArray::create<double4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneWeights");
vsiteLocalCoordsParams = CudaArray::create<double>(context, max(1, 12*numLocalCoords), "vsiteLocalCoordinatesParams");
vsiteLocalCoordsWeights = CudaArray::create<double>(context, max(1, (int) vsiteLocalCoordsWeightVec.size()), "vsiteLocalCoordsWeights");
vsiteLocalCoordsPos = CudaArray::create<double4>(context, max(1, (int) vsiteLocalCoordsPosVec.size()), "vsiteLocalCoordsPos");
if (num2Avg > 0)
vsite2AvgWeights->upload(vsite2AvgWeightVec);
if (num3Avg > 0)
vsite3AvgWeights->upload(vsite3AvgWeightVec);
if (numOutOfPlane > 0)
vsiteOutOfPlaneWeights->upload(vsiteOutOfPlaneWeightVec);
if (numLocalCoords > 0)
vsiteLocalCoordsParams->upload(vsiteLocalCoordsParamVec);
if (numLocalCoords > 0) {
vsiteLocalCoordsWeights->upload(vsiteLocalCoordsWeightVec);
vsiteLocalCoordsPos->upload(vsiteLocalCoordsPosVec);
}
}
else {
vsite2AvgWeights = CudaArray::create<float2>(context, max(1, num2Avg), "vsite2AvgWeights");
vsite3AvgWeights = CudaArray::create<float4>(context, max(1, num3Avg), "vsite3AvgWeights");
vsiteOutOfPlaneWeights = CudaArray::create<float4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneWeights");
vsiteLocalCoordsParams = CudaArray::create<float>(context, max(1, 12*numLocalCoords), "vsiteLocalCoordinatesParams");
vsiteLocalCoordsWeights = CudaArray::create<float>(context, max(1, (int) vsiteLocalCoordsWeightVec.size()), "vsiteLocalCoordsWeights");
vsiteLocalCoordsPos = CudaArray::create<float4>(context, max(1, (int) vsiteLocalCoordsPosVec.size()), "vsiteLocalCoordsPos");
if (num2Avg > 0) {
vector<float2> floatWeights(num2Avg);
for (int i = 0; i < num2Avg; i++)
......@@ -557,10 +569,14 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
vsiteOutOfPlaneWeights->upload(floatWeights);
}
if (numLocalCoords > 0) {
vector<float> floatParams(vsiteLocalCoordsParamVec.size());
for (int i = 0; i < (int) vsiteLocalCoordsParamVec.size(); i++)
floatParams[i] = (float) vsiteLocalCoordsParamVec[i];
vsiteLocalCoordsParams->upload(floatParams);
vector<float> floatWeights(vsiteLocalCoordsWeightVec.size());
for (int i = 0; i < (int) vsiteLocalCoordsWeightVec.size(); i++)
floatWeights[i] = (float) vsiteLocalCoordsWeightVec[i];
vsiteLocalCoordsWeights->upload(floatWeights);
vector<float4> floatPos(vsiteLocalCoordsPosVec.size());
for (int i = 0; i < (int) vsiteLocalCoordsPosVec.size(); i++)
floatPos[i] = make_float4((float) vsiteLocalCoordsPosVec[i].x, (float) vsiteLocalCoordsPosVec[i].y, (float) vsiteLocalCoordsPosVec[i].z, 0.0f);
vsiteLocalCoordsPos->upload(floatPos);
}
}
......@@ -644,10 +660,16 @@ CudaIntegrationUtilities::~CudaIntegrationUtilities() {
delete vsiteOutOfPlaneAtoms;
if (vsiteOutOfPlaneWeights != NULL)
delete vsiteOutOfPlaneWeights;
if (vsiteLocalCoordsIndex != NULL)
delete vsiteLocalCoordsIndex;
if (vsiteLocalCoordsAtoms != NULL)
delete vsiteLocalCoordsAtoms;
if (vsiteLocalCoordsParams != NULL)
delete vsiteLocalCoordsParams;
if (vsiteLocalCoordsWeights != NULL)
delete vsiteLocalCoordsWeights;
if (vsiteLocalCoordsPos != NULL)
delete vsiteLocalCoordsPos;
if (vsiteLocalCoordsStartIndex != NULL)
delete vsiteLocalCoordsStartIndex;
}
void CudaIntegrationUtilities::setNextStepSize(double size) {
......@@ -747,7 +769,9 @@ void CudaIntegrationUtilities::computeVirtualSites() {
void* args[] = {&context.getPosq().getDevicePointer(), &posCorrection, &vsite2AvgAtoms->getDevicePointer(), &vsite2AvgWeights->getDevicePointer(),
&vsite3AvgAtoms->getDevicePointer(), &vsite3AvgWeights->getDevicePointer(),
&vsiteOutOfPlaneAtoms->getDevicePointer(), &vsiteOutOfPlaneWeights->getDevicePointer(),
&vsiteLocalCoordsAtoms->getDevicePointer(), &vsiteLocalCoordsParams->getDevicePointer()};
&vsiteLocalCoordsIndex->getDevicePointer(), &vsiteLocalCoordsAtoms->getDevicePointer(),
&vsiteLocalCoordsWeights->getDevicePointer(), &vsiteLocalCoordsPos->getDevicePointer(),
&vsiteLocalCoordsStartIndex->getDevicePointer()};
context.executeKernel(vsitePositionKernel, args, numVsites);
}
}
......@@ -759,7 +783,9 @@ void CudaIntegrationUtilities::distributeForcesFromVirtualSites() {
&vsite2AvgAtoms->getDevicePointer(), &vsite2AvgWeights->getDevicePointer(),
&vsite3AvgAtoms->getDevicePointer(), &vsite3AvgWeights->getDevicePointer(),
&vsiteOutOfPlaneAtoms->getDevicePointer(), &vsiteOutOfPlaneWeights->getDevicePointer(),
&vsiteLocalCoordsAtoms->getDevicePointer(), &vsiteLocalCoordsParams->getDevicePointer()};
&vsiteLocalCoordsIndex->getDevicePointer(), &vsiteLocalCoordsAtoms->getDevicePointer(),
&vsiteLocalCoordsWeights->getDevicePointer(), &vsiteLocalCoordsPos->getDevicePointer(),
&vsiteLocalCoordsStartIndex->getDevicePointer()};
context.executeKernel(vsiteForceKernel, args, numVsites);
}
}
......
platforms/cuda/src/kernels/integrationUtilities.cu
View file @ 0827ad50
......@@ -680,7 +680,9 @@ extern "C" __global__ void updateCCMAAtomPositions(const int* __restrict__ numAt
extern "C" __global__ void computeVirtualSites(real4* __restrict__ posq, real4* __restrict__ posqCorrection, const int4* __restrict__ avg2Atoms, const real2* __restrict__ avg2Weights,
const int4* __restrict__ avg3Atoms, const real4* __restrict__ avg3Weights,
const int4* __restrict__ outOfPlaneAtoms, const real4* __restrict__ outOfPlaneWeights,
const int4* __restrict__ localCoordsAtoms, const real* __restrict__ localCoordsParams) {
const int* __restrict__ localCoordsIndex, const int* __restrict__ localCoordsAtoms,
const real* __restrict__ localCoordsWeights, const real4* __restrict__ localCoordsPos,
const int* __restrict__ localCoordsStartIndex) {
// Two particle average sites.
......@@ -732,30 +734,31 @@ extern "C" __global__ void computeVirtualSites(real4* __restrict__ posq, real4*
// Local coordinates sites.
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_LOCAL_COORDS; index += blockDim.x*gridDim.x) {
int4 atoms = localCoordsAtoms[index];
const real* params = &localCoordsParams[12*index];
mixed4 pos = loadPos(posq, posqCorrection, atoms.x);
mixed4 pos1_4 = loadPos(posq, posqCorrection, atoms.y);
mixed4 pos2_4 = loadPos(posq, posqCorrection, atoms.z);
mixed4 pos3_4 = loadPos(posq, posqCorrection, atoms.w);
mixed3 pos1 = make_mixed3(pos1_4.x, pos1_4.y, pos1_4.z);
mixed3 pos2 = make_mixed3(pos2_4.x, pos2_4.y, pos2_4.z);
mixed3 pos3 = make_mixed3(pos3_4.x, pos3_4.y, pos3_4.z);
mixed3 originWeights = make_mixed3(params[0], params[1], params[2]);
mixed3 xWeights = make_mixed3(params[3], params[4], params[5]);
mixed3 yWeights = make_mixed3(params[6], params[7], params[8]);
mixed3 localPosition = make_mixed3(params[9], params[10], params[11]);
mixed3 origin = pos1*originWeights.x + pos2*originWeights.y + pos3*originWeights.z;
mixed3 xdir = pos1*xWeights.x + pos2*xWeights.y + pos3*xWeights.z;
mixed3 ydir = pos1*yWeights.x + pos2*yWeights.y + pos3*yWeights.z;
int siteAtomIndex = localCoordsIndex[index];
int start = localCoordsStartIndex[index];
int end = localCoordsStartIndex[index+1];
mixed3 origin = make_mixed3(0), xdir = make_mixed3(0), ydir = make_mixed3(0);
for (int j = start; j < end; j++) {
mixed3 pos = trimTo3(loadPos(posq, posqCorrection, localCoordsAtoms[j]));
origin += pos*localCoordsWeights[3*j];
xdir += pos*localCoordsWeights[3*j+1];
ydir += pos*localCoordsWeights[3*j+2];
}
mixed3 zdir = cross(xdir, ydir);
xdir *= rsqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
zdir *= rsqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
mixed normXdir = sqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
mixed normZdir = sqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
mixed invNormXdir = (normXdir > 0 ? 1/normXdir : 0);
mixed invNormZdir = (normZdir > 0 ? 1/normZdir : 0);
xdir *= invNormXdir;
zdir *= invNormZdir;
ydir = cross(zdir, xdir);
real4 localPosition_4 = localCoordsPos[index];
mixed3 localPosition = make_mixed3(localPosition_4.x, localPosition_4.y, localPosition_4.z);
mixed4 pos = loadPos(posq, posqCorrection, siteAtomIndex);
pos.x = origin.x + xdir.x*localPosition.x + ydir.x*localPosition.y + zdir.x*localPosition.z;
pos.y = origin.y + xdir.y*localPosition.x + ydir.y*localPosition.y + zdir.y*localPosition.z;
pos.z = origin.z + xdir.z*localPosition.x + ydir.z*localPosition.y + zdir.z*localPosition.z;
storePos(posq, posqCorrection, atoms.x, pos);
storePos(posq, posqCorrection, siteAtomIndex, pos);
}
}
......@@ -778,7 +781,9 @@ extern "C" __global__ void distributeVirtualSiteForces(const real4* __restrict__
const int4* __restrict__ avg2Atoms, const real2* __restrict__ avg2Weights,
const int4* __restrict__ avg3Atoms, const real4* __restrict__ avg3Weights,
const int4* __restrict__ outOfPlaneAtoms, const real4* __restrict__ outOfPlaneWeights,
const int4* __restrict__ localCoordsAtoms, const real* __restrict__ localCoordsParams) {
const int* __restrict__ localCoordsIndex, const int* __restrict__ localCoordsAtoms,
const real* __restrict__ localCoordsWeights, const real4* __restrict__ localCoordsPos,
const int* __restrict__ localCoordsStartIndex) {
// Two particle average sites.
......@@ -826,87 +831,56 @@ extern "C" __global__ void distributeVirtualSiteForces(const real4* __restrict__
// Local coordinates sites.
for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_LOCAL_COORDS; index += blockDim.x*gridDim.x) {
int4 atoms = localCoordsAtoms[index];
const real* params = &localCoordsParams[12*index];
mixed4 pos = loadPos(posq, posqCorrection, atoms.x);
mixed4 pos1_4 = loadPos(posq, posqCorrection, atoms.y);
mixed4 pos2_4 = loadPos(posq, posqCorrection, atoms.z);
mixed4 pos3_4 = loadPos(posq, posqCorrection, atoms.w);
mixed3 pos1 = make_mixed3(pos1_4.x, pos1_4.y, pos1_4.z);
mixed3 pos2 = make_mixed3(pos2_4.x, pos2_4.y, pos2_4.z);
mixed3 pos3 = make_mixed3(pos3_4.x, pos3_4.y, pos3_4.z);
mixed3 originWeights = make_mixed3(params[0], params[1], params[2]);
mixed3 wx = make_mixed3(params[3], params[4], params[5]);
mixed3 wy = make_mixed3(params[6], params[7], params[8]);
mixed3 localPosition = make_mixed3(params[9], params[10], params[11]);
mixed3 origin = pos1*originWeights.x + pos2*originWeights.y + pos3*originWeights.z;
mixed3 xdir = pos1*wx.x + pos2*wx.y + pos3*wx.z;
mixed3 ydir = pos1*wy.x + pos2*wy.y + pos3*wy.z;
int siteAtomIndex = localCoordsIndex[index];
int start = localCoordsStartIndex[index];
int end = localCoordsStartIndex[index+1];
mixed3 origin = make_mixed3(0), xdir = make_mixed3(0), ydir = make_mixed3(0);
for (int j = start; j < end; j++) {
mixed3 pos = trimTo3(loadPos(posq, posqCorrection, localCoordsAtoms[j]));
origin += pos*localCoordsWeights[3*j];
xdir += pos*localCoordsWeights[3*j+1];
ydir += pos*localCoordsWeights[3*j+2];
}
mixed3 zdir = cross(xdir, ydir);
mixed invNormXdir = rsqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
mixed invNormZdir = rsqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
mixed normXdir = sqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
mixed normZdir = sqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
mixed invNormXdir = (normXdir > 0 ? 1/normXdir : 0);
mixed invNormZdir = (normZdir > 0 ? 1/normZdir : 0);
mixed3 dx = xdir*invNormXdir;
mixed3 dz = zdir*invNormZdir;
mixed3 dy = cross(dz, dx);
real4 localPosition_4 = localCoordsPos[index];
mixed3 localPosition = make_mixed3(localPosition_4.x, localPosition_4.y, localPosition_4.z);
// The derivatives for this case are very complicated. They were computed with SymPy then simplified by hand.
mixed t11 = (wx.x*ydir.x-wy.x*xdir.x)*invNormZdir;
mixed t12 = (wx.x*ydir.y-wy.x*xdir.y)*invNormZdir;
mixed t13 = (wx.x*ydir.z-wy.x*xdir.z)*invNormZdir;
mixed t21 = (wx.y*ydir.x-wy.y*xdir.x)*invNormZdir;
mixed t22 = (wx.y*ydir.y-wy.y*xdir.y)*invNormZdir;
mixed t23 = (wx.y*ydir.z-wy.y*xdir.z)*invNormZdir;
mixed t31 = (wx.z*ydir.x-wy.z*xdir.x)*invNormZdir;
mixed t32 = (wx.z*ydir.y-wy.z*xdir.y)*invNormZdir;
mixed t33 = (wx.z*ydir.z-wy.z*xdir.z)*invNormZdir;
mixed sx1 = t13*dz.y-t12*dz.z;
mixed sy1 = t11*dz.z-t13*dz.x;
mixed sz1 = t12*dz.x-t11*dz.y;
mixed sx2 = t23*dz.y-t22*dz.z;
mixed sy2 = t21*dz.z-t23*dz.x;
mixed sz2 = t22*dz.x-t21*dz.y;
mixed sx3 = t33*dz.y-t32*dz.z;
mixed sy3 = t31*dz.z-t33*dz.x;
mixed sz3 = t32*dz.x-t31*dz.y;
mixed3 wxScaled = wx*invNormXdir;
real3 f = loadForce(atoms.x, force);
real3 f = loadForce(siteAtomIndex, force);
mixed3 fp1 = localPosition*f.x;
mixed3 fp2 = localPosition*f.y;
mixed3 fp3 = localPosition*f.z;
real3 f1 = make_real3(0);
real3 f2 = make_real3(0);
real3 f3 = make_real3(0);
f1.x += fp1.x*wxScaled.x*(1-dx.x*dx.x) + fp1.z*(dz.x*sx1 ) + fp1.y*((-dx.x*dy.x )*wxScaled.x + dy.x*sx1 - dx.y*t12 - dx.z*t13) + f.x*originWeights.x;
f1.y += fp1.x*wxScaled.x*( -dx.x*dx.y) + fp1.z*(dz.x*sy1+t13) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.x + dy.x*sy1 + dx.y*t11);
f1.z += fp1.x*wxScaled.x*( -dx.x*dx.z) + fp1.z*(dz.x*sz1-t12) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.x + dy.x*sz1 + dx.z*t11);
f2.x += fp1.x*wxScaled.y*(1-dx.x*dx.x) + fp1.z*(dz.x*sx2 ) + fp1.y*((-dx.x*dy.x )*wxScaled.y + dy.x*sx2 - dx.y*t22 - dx.z*t23) + f.x*originWeights.y;
f2.y += fp1.x*wxScaled.y*( -dx.x*dx.y) + fp1.z*(dz.x*sy2+t23) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.y + dy.x*sy2 + dx.y*t21);
f2.z += fp1.x*wxScaled.y*( -dx.x*dx.z) + fp1.z*(dz.x*sz2-t22) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.y + dy.x*sz2 + dx.z*t21);
f3.x += fp1.x*wxScaled.z*(1-dx.x*dx.x) + fp1.z*(dz.x*sx3 ) + fp1.y*((-dx.x*dy.x )*wxScaled.z + dy.x*sx3 - dx.y*t32 - dx.z*t33) + f.x*originWeights.z;
f3.y += fp1.x*wxScaled.z*( -dx.x*dx.y) + fp1.z*(dz.x*sy3+t33) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.z + dy.x*sy3 + dx.y*t31);
f3.z += fp1.x*wxScaled.z*( -dx.x*dx.z) + fp1.z*(dz.x*sz3-t32) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.z + dy.x*sz3 + dx.z*t31);
f1.x += fp2.x*wxScaled.x*( -dx.y*dx.x) + fp2.z*(dz.y*sx1-t13) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.x - dy.y*sx1 - dx.x*t12);
f1.y += fp2.x*wxScaled.x*(1-dx.y*dx.y) + fp2.z*(dz.y*sy1 ) - fp2.y*(( dx.y*dy.y )*wxScaled.x - dy.y*sy1 + dx.x*t11 + dx.z*t13) + f.y*originWeights.x;
f1.z += fp2.x*wxScaled.x*( -dx.y*dx.z) + fp2.z*(dz.y*sz1+t11) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.x - dy.y*sz1 - dx.z*t12);
f2.x += fp2.x*wxScaled.y*( -dx.y*dx.x) + fp2.z*(dz.y*sx2-t23) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.y - dy.y*sx2 - dx.x*t22);
f2.y += fp2.x*wxScaled.y*(1-dx.y*dx.y) + fp2.z*(dz.y*sy2 ) - fp2.y*(( dx.y*dy.y )*wxScaled.y - dy.y*sy2 + dx.x*t21 + dx.z*t23) + f.y*originWeights.y;
f2.z += fp2.x*wxScaled.y*( -dx.y*dx.z) + fp2.z*(dz.y*sz2+t21) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.y - dy.y*sz2 - dx.z*t22);
f3.x += fp2.x*wxScaled.z*( -dx.y*dx.x) + fp2.z*(dz.y*sx3-t33) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.z - dy.y*sx3 - dx.x*t32);
f3.y += fp2.x*wxScaled.z*(1-dx.y*dx.y) + fp2.z*(dz.y*sy3 ) - fp2.y*(( dx.y*dy.y )*wxScaled.z - dy.y*sy3 + dx.x*t31 + dx.z*t33) + f.y*originWeights.z;
f3.z += fp2.x*wxScaled.z*( -dx.y*dx.z) + fp2.z*(dz.y*sz3+t31) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.z - dy.y*sz3 - dx.z*t32);
f1.x += fp3.x*wxScaled.x*( -dx.z*dx.x) + fp3.z*(dz.z*sx1+t12) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.x + dy.z*sx1 + dx.x*t13);
f1.y += fp3.x*wxScaled.x*( -dx.z*dx.y) + fp3.z*(dz.z*sy1-t11) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.x + dy.z*sy1 + dx.y*t13);
f1.z += fp3.x*wxScaled.x*(1-dx.z*dx.z) + fp3.z*(dz.z*sz1 ) + fp3.y*((-dx.z*dy.z )*wxScaled.x + dy.z*sz1 - dx.x*t11 - dx.y*t12) + f.z*originWeights.x;
f2.x += fp3.x*wxScaled.y*( -dx.z*dx.x) + fp3.z*(dz.z*sx2+t22) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.y + dy.z*sx2 + dx.x*t23);
f2.y += fp3.x*wxScaled.y*( -dx.z*dx.y) + fp3.z*(dz.z*sy2-t21) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.y + dy.z*sy2 + dx.y*t23);
f2.z += fp3.x*wxScaled.y*(1-dx.z*dx.z) + fp3.z*(dz.z*sz2 ) + fp3.y*((-dx.z*dy.z )*wxScaled.y + dy.z*sz2 - dx.x*t21 - dx.y*t22) + f.z*originWeights.y;
f3.x += fp3.x*wxScaled.z*( -dx.z*dx.x) + fp3.z*(dz.z*sx3+t32) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.z + dy.z*sx3 + dx.x*t33);
f3.y += fp3.x*wxScaled.z*( -dx.z*dx.y) + fp3.z*(dz.z*sy3-t31) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.z + dy.z*sy3 + dx.y*t33);
f3.z += fp3.x*wxScaled.z*(1-dx.z*dx.z) + fp3.z*(dz.z*sz3 ) + fp3.y*((-dx.z*dy.z )*wxScaled.z + dy.z*sz3 - dx.x*t31 - dx.y*t32) + f.z*originWeights.z;
addForce(atoms.y, force, f1);
addForce(atoms.z, force, f2);
addForce(atoms.w, force, f3);
for (int j = start; j < end; j++) {
real originWeight = localCoordsWeights[3*j];
real wx = localCoordsWeights[3*j+1];
real wy = localCoordsWeights[3*j+2];
mixed wxScaled = wx*invNormXdir;
mixed t1 = (wx*ydir.x-wy*xdir.x)*invNormZdir;
mixed t2 = (wx*ydir.y-wy*xdir.y)*invNormZdir;
mixed t3 = (wx*ydir.z-wy*xdir.z)*invNormZdir;
mixed sx = t3*dz.y-t2*dz.z;
mixed sy = t1*dz.z-t3*dz.x;
mixed sz = t2*dz.x-t1*dz.y;
real3 fresult = make_real3(0);
fresult.x += fp1.x*wxScaled*(1-dx.x*dx.x) + fp1.z*(dz.x*sx ) + fp1.y*((-dx.x*dy.x )*wxScaled + dy.x*sx - dx.y*t2 - dx.z*t3) + f.x*originWeight;
fresult.y += fp1.x*wxScaled*( -dx.x*dx.y) + fp1.z*(dz.x*sy+t3) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled + dy.x*sy + dx.y*t1);
fresult.z += fp1.x*wxScaled*( -dx.x*dx.z) + fp1.z*(dz.x*sz-t2) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled + dy.x*sz + dx.z*t1);
fresult.x += fp2.x*wxScaled*( -dx.y*dx.x) + fp2.z*(dz.y*sx-t3) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled - dy.y*sx - dx.x*t2);
fresult.y += fp2.x*wxScaled*(1-dx.y*dx.y) + fp2.z*(dz.y*sy ) - fp2.y*(( dx.y*dy.y )*wxScaled - dy.y*sy + dx.x*t1 + dx.z*t3) + f.y*originWeight;
fresult.z += fp2.x*wxScaled*( -dx.y*dx.z) + fp2.z*(dz.y*sz+t1) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled - dy.y*sz - dx.z*t2);
fresult.x += fp3.x*wxScaled*( -dx.z*dx.x) + fp3.z*(dz.z*sx+t2) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled + dy.z*sx + dx.x*t3);
fresult.y += fp3.x*wxScaled*( -dx.z*dx.y) + fp3.z*(dz.z*sy-t1) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled + dy.z*sy + dx.y*t3);
fresult.z += fp3.x*wxScaled*(1-dx.z*dx.z) + fp3.z*(dz.z*sz ) + fp3.y*((-dx.z*dy.z )*wxScaled + dy.z*sz - dx.x*t1 - dx.y*t2) + f.z*originWeight;
addForce(localCoordsAtoms[j], force, fresult);
}
}
}
......
platforms/opencl/include/OpenCLIntegrationUtilities.h
View file @ 0827ad50
......@@ -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) 2009-2014 Stanford University and the Authors. *
* Portions copyright (c) 2009-2017 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -156,8 +156,11 @@ private:
OpenCLArray* vsite3AvgWeights;
OpenCLArray* vsiteOutOfPlaneAtoms;
OpenCLArray* vsiteOutOfPlaneWeights;
OpenCLArray* vsiteLocalCoordsIndex;
OpenCLArray* vsiteLocalCoordsAtoms;
OpenCLArray* vsiteLocalCoordsParams;
OpenCLArray* vsiteLocalCoordsWeights;
OpenCLArray* vsiteLocalCoordsPos;
OpenCLArray* vsiteLocalCoordsStartIndex;
int randomPos;
int lastSeed, numVsites;
bool hasInitializedPosConstraintKernels, hasInitializedVelConstraintKernels, ccmaUseDirectBuffer, hasOverlappingVsites;
......
platforms/opencl/src/OpenCLIntegrationUtilities.cpp
View file @ 0827ad50
......@@ -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) 2009-2015 Stanford University and the Authors. *
* Portions copyright (c) 2009-2017 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -102,7 +102,8 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
ccmaReducedMass(NULL), ccmaAtomConstraints(NULL), ccmaNumAtomConstraints(NULL), ccmaConstraintMatrixColumn(NULL),
ccmaConstraintMatrixValue(NULL), ccmaDelta1(NULL), ccmaDelta2(NULL), ccmaConverged(NULL), ccmaConvergedHostBuffer(NULL),
vsite2AvgAtoms(NULL), vsite2AvgWeights(NULL), vsite3AvgAtoms(NULL), vsite3AvgWeights(NULL),
vsiteOutOfPlaneAtoms(NULL), vsiteOutOfPlaneWeights(NULL), vsiteLocalCoordsAtoms(NULL), vsiteLocalCoordsParams(NULL),
vsiteOutOfPlaneAtoms(NULL), vsiteOutOfPlaneWeights(NULL), vsiteLocalCoordsIndex(NULL), vsiteLocalCoordsAtoms(NULL),
vsiteLocalCoordsWeights(NULL), vsiteLocalCoordsPos(NULL), vsiteLocalCoordsStartIndex(NULL),
hasInitializedPosConstraintKernels(false), hasInitializedVelConstraintKernels(false), hasOverlappingVsites(false) {
// Create workspace arrays.
......@@ -497,8 +498,11 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
vector<mm_double4> vsite3AvgWeightVec;
vector<mm_int4> vsiteOutOfPlaneAtomVec;
vector<mm_double4> vsiteOutOfPlaneWeightVec;
vector<mm_int4> vsiteLocalCoordsAtomVec;
vector<cl_double> vsiteLocalCoordsParamVec;
vector<cl_int> vsiteLocalCoordsIndexVec;
vector<cl_int> vsiteLocalCoordsAtomVec;
vector<cl_int> vsiteLocalCoordsStartVec;
vector<cl_double> vsiteLocalCoordsWeightVec;
vector<mm_double4> vsiteLocalCoordsPosVec;
for (int i = 0; i < numAtoms; i++) {
if (system.isVirtualSite(i)) {
if (dynamic_cast<const TwoParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
......@@ -523,65 +527,73 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
vsiteOutOfPlaneWeightVec.push_back(mm_double4(site.getWeight12(), site.getWeight13(), site.getWeightCross(), 0.0));
}
else if (dynamic_cast<const LocalCoordinatesSite*>(&system.getVirtualSite(i)) != NULL) {
// An out of plane site.
// A local coordinates site.
const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
vsiteLocalCoordsAtomVec.push_back(mm_int4(i, site.getParticle(0), site.getParticle(1), site.getParticle(2)));
Vec3 origin = site.getOriginWeights();
Vec3 x = site.getXWeights();
Vec3 y = site.getYWeights();
int numParticles = site.getNumParticles();
vector<double> origin, x, y;
site.getOriginWeights(origin);
site.getXWeights(x);
site.getYWeights(y);
vsiteLocalCoordsIndexVec.push_back(i);
vsiteLocalCoordsStartVec.push_back(vsiteLocalCoordsAtomVec.size());
for (int j = 0; j < numParticles; j++) {
vsiteLocalCoordsAtomVec.push_back(site.getParticle(j));
vsiteLocalCoordsWeightVec.push_back(origin[j]);
vsiteLocalCoordsWeightVec.push_back(x[j]);
vsiteLocalCoordsWeightVec.push_back(y[j]);
}
Vec3 pos = site.getLocalPosition();
vsiteLocalCoordsParamVec.push_back(origin[0]);
vsiteLocalCoordsParamVec.push_back(origin[1]);
vsiteLocalCoordsParamVec.push_back(origin[2]);
vsiteLocalCoordsParamVec.push_back(x[0]);
vsiteLocalCoordsParamVec.push_back(x[1]);
vsiteLocalCoordsParamVec.push_back(x[2]);
vsiteLocalCoordsParamVec.push_back(y[0]);
vsiteLocalCoordsParamVec.push_back(y[1]);
vsiteLocalCoordsParamVec.push_back(y[2]);
vsiteLocalCoordsParamVec.push_back(pos[0]);
vsiteLocalCoordsParamVec.push_back(pos[1]);
vsiteLocalCoordsParamVec.push_back(pos[2]);
vsiteLocalCoordsPosVec.push_back(mm_double4(pos[0], pos[1], pos[2], 0.0));
}
}
}
vsiteLocalCoordsStartVec.push_back(vsiteLocalCoordsAtomVec.size());
int num2Avg = vsite2AvgAtomVec.size();
int num3Avg = vsite3AvgAtomVec.size();
int numOutOfPlane = vsiteOutOfPlaneAtomVec.size();
int numLocalCoords = vsiteLocalCoordsAtomVec.size();
int numLocalCoords = vsiteLocalCoordsPosVec.size();
numVsites = num2Avg+num3Avg+numOutOfPlane+numLocalCoords;
vsite2AvgAtoms = OpenCLArray::create<mm_int4>(context, max(1, num2Avg), "vsite2AvgAtoms");
vsite3AvgAtoms = OpenCLArray::create<mm_int4>(context, max(1, num3Avg), "vsite3AvgAtoms");
vsiteOutOfPlaneAtoms = OpenCLArray::create<mm_int4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneAtoms");
vsiteLocalCoordsAtoms = OpenCLArray::create<mm_int4>(context, max(1, numLocalCoords), "vsiteLocalCoordinatesAtoms");
vsiteLocalCoordsIndex = OpenCLArray::create<cl_int>(context, max(1, (int) vsiteLocalCoordsIndexVec.size()), "vsiteLocalCoordsIndex");
vsiteLocalCoordsAtoms = OpenCLArray::create<cl_int>(context, max(1, (int) vsiteLocalCoordsAtomVec.size()), "vsiteLocalCoordsAtoms");
vsiteLocalCoordsStartIndex = OpenCLArray::create<cl_int>(context, max(1, (int) vsiteLocalCoordsStartVec.size()), "vsiteLocalCoordsStartIndex");
if (num2Avg > 0)
vsite2AvgAtoms->upload(vsite2AvgAtomVec);
if (num3Avg > 0)
vsite3AvgAtoms->upload(vsite3AvgAtomVec);
if (numOutOfPlane > 0)
vsiteOutOfPlaneAtoms->upload(vsiteOutOfPlaneAtomVec);
if (numLocalCoords > 0)
if (numLocalCoords > 0) {
vsiteLocalCoordsIndex->upload(vsiteLocalCoordsIndexVec);
vsiteLocalCoordsAtoms->upload(vsiteLocalCoordsAtomVec);
vsiteLocalCoordsStartIndex->upload(vsiteLocalCoordsStartVec);
}
if (context.getUseDoublePrecision()) {
vsite2AvgWeights = OpenCLArray::create<mm_double2>(context, max(1, num2Avg), "vsite2AvgWeights");
vsite3AvgWeights = OpenCLArray::create<mm_double4>(context, max(1, num3Avg), "vsite3AvgWeights");
vsiteOutOfPlaneWeights = OpenCLArray::create<mm_double4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneWeights");
vsiteLocalCoordsParams = OpenCLArray::create<cl_double>(context, max(1, 12*numLocalCoords), "vsiteLocalCoordinatesParams");
vsiteLocalCoordsWeights = OpenCLArray::create<cl_double>(context, max(1, (int) vsiteLocalCoordsWeightVec.size()), "vsiteLocalCoordsWeights");
vsiteLocalCoordsPos = OpenCLArray::create<mm_double4>(context, max(1, (int) vsiteLocalCoordsPosVec.size()), "vsiteLocalCoordsPos");
if (num2Avg > 0)
vsite2AvgWeights->upload(vsite2AvgWeightVec);
if (num3Avg > 0)
vsite3AvgWeights->upload(vsite3AvgWeightVec);
if (numOutOfPlane > 0)
vsiteOutOfPlaneWeights->upload(vsiteOutOfPlaneWeightVec);
if (numLocalCoords > 0)
vsiteLocalCoordsParams->upload(vsiteLocalCoordsParamVec);
if (numLocalCoords > 0) {
vsiteLocalCoordsWeights->upload(vsiteLocalCoordsWeightVec);
vsiteLocalCoordsPos->upload(vsiteLocalCoordsPosVec);
}
}
else {
vsite2AvgWeights = OpenCLArray::create<mm_float2>(context, max(1, num2Avg), "vsite2AvgWeights");
vsite3AvgWeights = OpenCLArray::create<mm_float4>(context, max(1, num3Avg), "vsite3AvgWeights");
vsiteOutOfPlaneWeights = OpenCLArray::create<mm_float4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneWeights");
vsiteLocalCoordsParams = OpenCLArray::create<float>(context, max(1, 12*numLocalCoords), "vsiteLocalCoordinatesParams");
vsiteLocalCoordsWeights = OpenCLArray::create<cl_float>(context, max(1, (int) vsiteLocalCoordsWeightVec.size()), "vsiteLocalCoordsWeights");
vsiteLocalCoordsPos = OpenCLArray::create<mm_float4>(context, max(1, (int) vsiteLocalCoordsPosVec.size()), "vsiteLocalCoordsPos");
if (num2Avg > 0) {
vector<mm_float2> floatWeights(num2Avg);
for (int i = 0; i < num2Avg; i++)
......@@ -601,10 +613,14 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
vsiteOutOfPlaneWeights->upload(floatWeights);
}
if (numLocalCoords > 0) {
vector<cl_float> floatParams(vsiteLocalCoordsParamVec.size());
for (int i = 0; i < (int) vsiteLocalCoordsParamVec.size(); i++)
floatParams[i] = (cl_float) vsiteLocalCoordsParamVec[i];
vsiteLocalCoordsParams->upload(floatParams);
vector<cl_float> floatWeights(vsiteLocalCoordsWeightVec.size());
for (int i = 0; i < (int) vsiteLocalCoordsWeightVec.size(); i++)
floatWeights[i] = (cl_float) vsiteLocalCoordsWeightVec[i];
vsiteLocalCoordsWeights->upload(floatWeights);
vector<mm_float4> floatPos(vsiteLocalCoordsPosVec.size());
for (int i = 0; i < (int) vsiteLocalCoordsPosVec.size(); i++)
floatPos[i] = mm_float4((float) vsiteLocalCoordsPosVec[i].x, (float) vsiteLocalCoordsPosVec[i].y, (float) vsiteLocalCoordsPosVec[i].z, 0.0f);
vsiteLocalCoordsPos->upload(floatPos);
}
}
......@@ -645,8 +661,11 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
vsitePositionKernel.setArg<cl::Buffer>(index++, vsite3AvgWeights->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteOutOfPlaneAtoms->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteOutOfPlaneWeights->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsIndex->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsAtoms->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsParams->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsWeights->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsPos->getDeviceBuffer());
vsitePositionKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsStartIndex->getDeviceBuffer());
vsiteForceKernel = cl::Kernel(vsiteProgram, "distributeForces");
index = 0;
vsiteForceKernel.setArg<cl::Buffer>(index++, context.getPosq().getDeviceBuffer());
......@@ -661,8 +680,11 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
vsiteForceKernel.setArg<cl::Buffer>(index++, vsite3AvgWeights->getDeviceBuffer());
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteOutOfPlaneAtoms->getDeviceBuffer());
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteOutOfPlaneWeights->getDeviceBuffer());
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsIndex->getDeviceBuffer());
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsAtoms->getDeviceBuffer());
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsParams->getDeviceBuffer());
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsWeights->getDeviceBuffer());
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsPos->getDeviceBuffer());
vsiteForceKernel.setArg<cl::Buffer>(index++, vsiteLocalCoordsStartIndex->getDeviceBuffer());
if (hasOverlappingVsites && context.getSupports64BitGlobalAtomics())
vsiteAddForcesKernel = cl::Kernel(vsiteProgram, "addDistributedForces");
}
......@@ -718,10 +740,16 @@ OpenCLIntegrationUtilities::~OpenCLIntegrationUtilities() {
delete vsiteOutOfPlaneAtoms;
if (vsiteOutOfPlaneWeights != NULL)
delete vsiteOutOfPlaneWeights;
if (vsiteLocalCoordsIndex != NULL)
delete vsiteLocalCoordsIndex;
if (vsiteLocalCoordsAtoms != NULL)
delete vsiteLocalCoordsAtoms;
if (vsiteLocalCoordsParams != NULL)
delete vsiteLocalCoordsParams;
if (vsiteLocalCoordsWeights != NULL)
delete vsiteLocalCoordsWeights;
if (vsiteLocalCoordsPos != NULL)
delete vsiteLocalCoordsPos;
if (vsiteLocalCoordsStartIndex != NULL)
delete vsiteLocalCoordsStartIndex;
}
void OpenCLIntegrationUtilities::setNextStepSize(double size) {
......
platforms/opencl/src/kernels/virtualSites.cl
View file @ 0827ad50
......@@ -33,7 +33,9 @@ __kernel void computeVirtualSites(__global real4* restrict posq,
__global const int4* restrict avg2Atoms, __global const real2* restrict avg2Weights,
__global const int4* restrict avg3Atoms, __global const real4* restrict avg3Weights,
__global const int4* restrict outOfPlaneAtoms, __global const real4* restrict outOfPlaneWeights,
__global const int4* restrict localCoordsAtoms, __global const real* restrict localCoordsParams) {
__global const int* restrict localCoordsIndex, __global const int* restrict localCoordsAtoms,
__global const real* restrict localCoordsWeights, __global const real4* restrict localCoordsPos,
__global const int* restrict localCoordsStartIndex) {
#ifndef USE_MIXED_PRECISION
__global real4* posqCorrection = 0;
#endif
......@@ -81,30 +83,30 @@ __kernel void computeVirtualSites(__global real4* restrict posq,
// Local coordinates sites.
for (int index = get_global_id(0); index < NUM_LOCAL_COORDS; index += get_global_size(0)) {
int4 atoms = localCoordsAtoms[index];
__global const real* params = &localCoordsParams[12*index];
mixed4 pos = loadPos(posq, posqCorrection, atoms.x);
mixed4 pos1_4 = loadPos(posq, posqCorrection, atoms.y);
mixed4 pos2_4 = loadPos(posq, posqCorrection, atoms.z);
mixed4 pos3_4 = loadPos(posq, posqCorrection, atoms.w);
mixed4 pos1 = (mixed4) (pos1_4.x, pos1_4.y, pos1_4.z, 0);
mixed4 pos2 = (mixed4) (pos2_4.x, pos2_4.y, pos2_4.z, 0);
mixed4 pos3 = (mixed4) (pos3_4.x, pos3_4.y, pos3_4.z, 0);
mixed4 originWeights = (mixed4) (params[0], params[1], params[2], 0);
mixed4 xWeights = (mixed4) (params[3], params[4], params[5], 0);
mixed4 yWeights = (mixed4) (params[6], params[7], params[8], 0);
mixed4 localPosition = (mixed4) (params[9], params[10], params[11], 0);
mixed4 origin = pos1*originWeights.x + pos2*originWeights.y + pos3*originWeights.z;
mixed4 xdir = pos1*xWeights.x + pos2*xWeights.y + pos3*xWeights.z;
mixed4 ydir = pos1*yWeights.x + pos2*yWeights.y + pos3*yWeights.z;
mixed4 zdir = cross(xdir, ydir);
xdir *= rsqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
zdir *= rsqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
int siteAtomIndex = localCoordsIndex[index];
int start = localCoordsStartIndex[index];
int end = localCoordsStartIndex[index+1];
mixed3 origin = 0, xdir = 0, ydir = 0;
for (int j = start; j < end; j++) {
mixed3 pos = loadPos(posq, posqCorrection, localCoordsAtoms[j]).xyz;
origin += pos*localCoordsWeights[3*j];
xdir += pos*localCoordsWeights[3*j+1];
ydir += pos*localCoordsWeights[3*j+2];
}
mixed3 zdir = cross(xdir, ydir);
mixed normXdir = sqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
mixed normZdir = sqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
mixed invNormXdir = (normXdir > 0 ? 1/normXdir : 0);
mixed invNormZdir = (normZdir > 0 ? 1/normZdir : 0);
xdir *= invNormXdir;
zdir *= invNormZdir;
ydir = cross(zdir, xdir);
mixed3 localPosition = convert_mixed4(localCoordsPos[index]).xyz;
mixed4 pos = loadPos(posq, posqCorrection, siteAtomIndex);
pos.x = origin.x + xdir.x*localPosition.x + ydir.x*localPosition.y + zdir.x*localPosition.z;
pos.y = origin.y + xdir.y*localPosition.x + ydir.y*localPosition.y + zdir.y*localPosition.z;
pos.z = origin.z + xdir.z*localPosition.x + ydir.z*localPosition.y + zdir.z*localPosition.z;
storePos(posq, posqCorrection, atoms.x, pos);
storePos(posq, posqCorrection, siteAtomIndex, pos);
}
}
......@@ -174,7 +176,9 @@ __kernel void distributeForces(__global const real4* restrict posq, __global rea
__global const int4* restrict avg2Atoms, __global const real2* restrict avg2Weights,
__global const int4* restrict avg3Atoms, __global const real4* restrict avg3Weights,
__global const int4* restrict outOfPlaneAtoms, __global const real4* restrict outOfPlaneWeights,
__global const int4* restrict localCoordsAtoms, __global const real* restrict localCoordsParams) {
__global const int* restrict localCoordsIndex, __global const int* restrict localCoordsAtoms,
__global const real* restrict localCoordsWeights, __global const real4* restrict localCoordsPos,
__global const int* restrict localCoordsStartIndex) {
#ifndef USE_MIXED_PRECISION
__global real4* posqCorrection = 0;
#endif
......@@ -225,86 +229,54 @@ __kernel void distributeForces(__global const real4* restrict posq, __global rea
// Local coordinates sites.
for (int index = get_global_id(0); index < NUM_LOCAL_COORDS; index += get_global_size(0)) {
int4 atoms = localCoordsAtoms[index];
__global const real* params = &localCoordsParams[12*index];
mixed4 pos = loadPos(posq, posqCorrection, atoms.x);
mixed4 pos1_4 = loadPos(posq, posqCorrection, atoms.y);
mixed4 pos2_4 = loadPos(posq, posqCorrection, atoms.z);
mixed4 pos3_4 = loadPos(posq, posqCorrection, atoms.w);
mixed4 pos1 = (mixed4) (pos1_4.x, pos1_4.y, pos1_4.z, 0);
mixed4 pos2 = (mixed4) (pos2_4.x, pos2_4.y, pos2_4.z, 0);
mixed4 pos3 = (mixed4) (pos3_4.x, pos3_4.y, pos3_4.z, 0);
mixed4 originWeights = (mixed4) (params[0], params[1], params[2], 0);
mixed4 wx = (mixed4) (params[3], params[4], params[5], 0);
mixed4 wy = (mixed4) (params[6], params[7], params[8], 0);
mixed4 localPosition = (mixed4) (params[9], params[10], params[11], 0);
mixed4 origin = pos1*originWeights.x + pos2*originWeights.y + pos3*originWeights.z;
mixed4 xdir = pos1*wx.x + pos2*wx.y + pos3*wx.z;
mixed4 ydir = pos1*wy.x + pos2*wy.y + pos3*wy.z;
mixed4 zdir = cross(xdir, ydir);
mixed invNormXdir = rsqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
mixed invNormZdir = rsqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
mixed4 dx = xdir*invNormXdir;
mixed4 dz = zdir*invNormZdir;
mixed4 dy = cross(dz, dx);
int siteAtomIndex = localCoordsIndex[index];
int start = localCoordsStartIndex[index];
int end = localCoordsStartIndex[index+1];
mixed3 origin = 0, xdir = 0, ydir = 0;
for (int j = start; j < end; j++) {
mixed3 pos = loadPos(posq, posqCorrection, localCoordsAtoms[j]).xyz;
origin += pos*localCoordsWeights[3*j];
xdir += pos*localCoordsWeights[3*j+1];
ydir += pos*localCoordsWeights[3*j+2];
}
mixed3 zdir = cross(xdir, ydir);
mixed normXdir = sqrt(xdir.x*xdir.x+xdir.y*xdir.y+xdir.z*xdir.z);
mixed normZdir = sqrt(zdir.x*zdir.x+zdir.y*zdir.y+zdir.z*zdir.z);
mixed invNormXdir = (normXdir > 0 ? 1/normXdir : 0);
mixed invNormZdir = (normZdir > 0 ? 1/normZdir : 0);
mixed3 dx = xdir*invNormXdir;
mixed3 dz = zdir*invNormZdir;
mixed3 dy = cross(dz, dx);
mixed3 localPosition = convert_mixed4(localCoordsPos[index]).xyz;
// The derivatives for this case are very complicated. They were computed with SymPy then simplified by hand.
mixed t11 = (wx.x*ydir.x-wy.x*xdir.x)*invNormZdir;
mixed t12 = (wx.x*ydir.y-wy.x*xdir.y)*invNormZdir;
mixed t13 = (wx.x*ydir.z-wy.x*xdir.z)*invNormZdir;
mixed t21 = (wx.y*ydir.x-wy.y*xdir.x)*invNormZdir;
mixed t22 = (wx.y*ydir.y-wy.y*xdir.y)*invNormZdir;
mixed t23 = (wx.y*ydir.z-wy.y*xdir.z)*invNormZdir;
mixed t31 = (wx.z*ydir.x-wy.z*xdir.x)*invNormZdir;
mixed t32 = (wx.z*ydir.y-wy.z*xdir.y)*invNormZdir;
mixed t33 = (wx.z*ydir.z-wy.z*xdir.z)*invNormZdir;
mixed sx1 = t13*dz.y-t12*dz.z;
mixed sy1 = t11*dz.z-t13*dz.x;
mixed sz1 = t12*dz.x-t11*dz.y;
mixed sx2 = t23*dz.y-t22*dz.z;
mixed sy2 = t21*dz.z-t23*dz.x;
mixed sz2 = t22*dz.x-t21*dz.y;
mixed sx3 = t33*dz.y-t32*dz.z;
mixed sy3 = t31*dz.z-t33*dz.x;
mixed sz3 = t32*dz.x-t31*dz.y;
mixed4 wxScaled = wx*invNormXdir;
real4 f = force[atoms.x];
real4 f1 = 0;
real4 f2 = 0;
real4 f3 = 0;
mixed4 fp1 = localPosition*f.x;
mixed4 fp2 = localPosition*f.y;
mixed4 fp3 = localPosition*f.z;
f1.x += fp1.x*wxScaled.x*(1-dx.x*dx.x) + fp1.z*(dz.x*sx1 ) + fp1.y*((-dx.x*dy.x )*wxScaled.x + dy.x*sx1 - dx.y*t12 - dx.z*t13) + f.x*originWeights.x;
f1.y += fp1.x*wxScaled.x*( -dx.x*dx.y) + fp1.z*(dz.x*sy1+t13) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.x + dy.x*sy1 + dx.y*t11);
f1.z += fp1.x*wxScaled.x*( -dx.x*dx.z) + fp1.z*(dz.x*sz1-t12) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.x + dy.x*sz1 + dx.z*t11);
f2.x += fp1.x*wxScaled.y*(1-dx.x*dx.x) + fp1.z*(dz.x*sx2 ) + fp1.y*((-dx.x*dy.x )*wxScaled.y + dy.x*sx2 - dx.y*t22 - dx.z*t23) + f.x*originWeights.y;
f2.y += fp1.x*wxScaled.y*( -dx.x*dx.y) + fp1.z*(dz.x*sy2+t23) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.y + dy.x*sy2 + dx.y*t21);
f2.z += fp1.x*wxScaled.y*( -dx.x*dx.z) + fp1.z*(dz.x*sz2-t22) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.y + dy.x*sz2 + dx.z*t21);
f3.x += fp1.x*wxScaled.z*(1-dx.x*dx.x) + fp1.z*(dz.x*sx3 ) + fp1.y*((-dx.x*dy.x )*wxScaled.z + dy.x*sx3 - dx.y*t32 - dx.z*t33) + f.x*originWeights.z;
f3.y += fp1.x*wxScaled.z*( -dx.x*dx.y) + fp1.z*(dz.x*sy3+t33) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled.z + dy.x*sy3 + dx.y*t31);
f3.z += fp1.x*wxScaled.z*( -dx.x*dx.z) + fp1.z*(dz.x*sz3-t32) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled.z + dy.x*sz3 + dx.z*t31);
f1.x += fp2.x*wxScaled.x*( -dx.y*dx.x) + fp2.z*(dz.y*sx1-t13) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.x - dy.y*sx1 - dx.x*t12);
f1.y += fp2.x*wxScaled.x*(1-dx.y*dx.y) + fp2.z*(dz.y*sy1 ) - fp2.y*(( dx.y*dy.y )*wxScaled.x - dy.y*sy1 + dx.x*t11 + dx.z*t13) + f.y*originWeights.x;
f1.z += fp2.x*wxScaled.x*( -dx.y*dx.z) + fp2.z*(dz.y*sz1+t11) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.x - dy.y*sz1 - dx.z*t12);
f2.x += fp2.x*wxScaled.y*( -dx.y*dx.x) + fp2.z*(dz.y*sx2-t23) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.y - dy.y*sx2 - dx.x*t22);
f2.y += fp2.x*wxScaled.y*(1-dx.y*dx.y) + fp2.z*(dz.y*sy2 ) - fp2.y*(( dx.y*dy.y )*wxScaled.y - dy.y*sy2 + dx.x*t21 + dx.z*t23) + f.y*originWeights.y;
f2.z += fp2.x*wxScaled.y*( -dx.y*dx.z) + fp2.z*(dz.y*sz2+t21) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.y - dy.y*sz2 - dx.z*t22);
f3.x += fp2.x*wxScaled.z*( -dx.y*dx.x) + fp2.z*(dz.y*sx3-t33) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled.z - dy.y*sx3 - dx.x*t32);
f3.y += fp2.x*wxScaled.z*(1-dx.y*dx.y) + fp2.z*(dz.y*sy3 ) - fp2.y*(( dx.y*dy.y )*wxScaled.z - dy.y*sy3 + dx.x*t31 + dx.z*t33) + f.y*originWeights.z;
f3.z += fp2.x*wxScaled.z*( -dx.y*dx.z) + fp2.z*(dz.y*sz3+t31) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled.z - dy.y*sz3 - dx.z*t32);
f1.x += fp3.x*wxScaled.x*( -dx.z*dx.x) + fp3.z*(dz.z*sx1+t12) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.x + dy.z*sx1 + dx.x*t13);
f1.y += fp3.x*wxScaled.x*( -dx.z*dx.y) + fp3.z*(dz.z*sy1-t11) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.x + dy.z*sy1 + dx.y*t13);
f1.z += fp3.x*wxScaled.x*(1-dx.z*dx.z) + fp3.z*(dz.z*sz1 ) + fp3.y*((-dx.z*dy.z )*wxScaled.x + dy.z*sz1 - dx.x*t11 - dx.y*t12) + f.z*originWeights.x;
f2.x += fp3.x*wxScaled.y*( -dx.z*dx.x) + fp3.z*(dz.z*sx2+t22) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.y + dy.z*sx2 + dx.x*t23);
f2.y += fp3.x*wxScaled.y*( -dx.z*dx.y) + fp3.z*(dz.z*sy2-t21) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.y + dy.z*sy2 + dx.y*t23);
f2.z += fp3.x*wxScaled.y*(1-dx.z*dx.z) + fp3.z*(dz.z*sz2 ) + fp3.y*((-dx.z*dy.z )*wxScaled.y + dy.z*sz2 - dx.x*t21 - dx.y*t22) + f.z*originWeights.y;
f3.x += fp3.x*wxScaled.z*( -dx.z*dx.x) + fp3.z*(dz.z*sx3+t32) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled.z + dy.z*sx3 + dx.x*t33);
f3.y += fp3.x*wxScaled.z*( -dx.z*dx.y) + fp3.z*(dz.z*sy3-t31) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled.z + dy.z*sy3 + dx.y*t33);
f3.z += fp3.x*wxScaled.z*(1-dx.z*dx.z) + fp3.z*(dz.z*sz3 ) + fp3.y*((-dx.z*dy.z )*wxScaled.z + dy.z*sz3 - dx.x*t31 - dx.y*t32) + f.z*originWeights.z;
ADD_FORCE(atoms.y, f1);
ADD_FORCE(atoms.z, f2);
ADD_FORCE(atoms.w, f3);
real4 f = force[siteAtomIndex];
mixed3 fp1 = localPosition*f.x;
mixed3 fp2 = localPosition*f.y;
mixed3 fp3 = localPosition*f.z;
for (int j = start; j < end; j++) {
real originWeight = localCoordsWeights[3*j];
real wx = localCoordsWeights[3*j+1];
real wy = localCoordsWeights[3*j+2];
mixed wxScaled = wx*invNormXdir;
mixed t1 = (wx*ydir.x-wy*xdir.x)*invNormZdir;
mixed t2 = (wx*ydir.y-wy*xdir.y)*invNormZdir;
mixed t3 = (wx*ydir.z-wy*xdir.z)*invNormZdir;
mixed sx = t3*dz.y-t2*dz.z;
mixed sy = t1*dz.z-t3*dz.x;
mixed sz = t2*dz.x-t1*dz.y;
real4 fresult = 0;
fresult.x += fp1.x*wxScaled*(1-dx.x*dx.x) + fp1.z*(dz.x*sx ) + fp1.y*((-dx.x*dy.x )*wxScaled + dy.x*sx - dx.y*t2 - dx.z*t3) + f.x*originWeight;
fresult.y += fp1.x*wxScaled*( -dx.x*dx.y) + fp1.z*(dz.x*sy+t3) + fp1.y*((-dx.y*dy.x-dz.z)*wxScaled + dy.x*sy + dx.y*t1);
fresult.z += fp1.x*wxScaled*( -dx.x*dx.z) + fp1.z*(dz.x*sz-t2) + fp1.y*((-dx.z*dy.x+dz.y)*wxScaled + dy.x*sz + dx.z*t1);
fresult.x += fp2.x*wxScaled*( -dx.y*dx.x) + fp2.z*(dz.y*sx-t3) - fp2.y*(( dx.x*dy.y-dz.z)*wxScaled - dy.y*sx - dx.x*t2);
fresult.y += fp2.x*wxScaled*(1-dx.y*dx.y) + fp2.z*(dz.y*sy ) - fp2.y*(( dx.y*dy.y )*wxScaled - dy.y*sy + dx.x*t1 + dx.z*t3) + f.y*originWeight;
fresult.z += fp2.x*wxScaled*( -dx.y*dx.z) + fp2.z*(dz.y*sz+t1) - fp2.y*(( dx.z*dy.y+dz.x)*wxScaled - dy.y*sz - dx.z*t2);
fresult.x += fp3.x*wxScaled*( -dx.z*dx.x) + fp3.z*(dz.z*sx+t2) + fp3.y*((-dx.x*dy.z-dz.y)*wxScaled + dy.z*sx + dx.x*t3);
fresult.y += fp3.x*wxScaled*( -dx.z*dx.y) + fp3.z*(dz.z*sy-t1) + fp3.y*((-dx.y*dy.z+dz.x)*wxScaled + dy.z*sy + dx.y*t3);
fresult.z += fp3.x*wxScaled*(1-dx.z*dx.z) + fp3.z*(dz.z*sz ) + fp3.y*((-dx.z*dy.z )*wxScaled + dy.z*sz - dx.x*t1 - dx.y*t2) + f.z*originWeight;
ADD_FORCE(localCoordsAtoms[j], fresult);
}
}
}
platforms/reference/src/SimTKReference/ReferenceVirtualSites.cpp
View file @ 0827ad50
......@@ -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) 2012-2014 Stanford University and the Authors. *
* Portions copyright (c) 2012-2017 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -70,22 +70,30 @@ void ReferenceVirtualSites::computePositions(const OpenMM::System& system, vecto
// A local coordinates site.
const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
Vec3 originWeights = site.getOriginWeights();
Vec3 xWeights = site.getXWeights();
Vec3 yWeights = site.getYWeights();
int numParticles = site.getNumParticles();
vector<double> originWeights, xWeights, yWeights;
site.getOriginWeights(originWeights);
site.getXWeights(xWeights);
site.getYWeights(yWeights);
Vec3 origin, xdir, ydir;
for (int j = 0; j < numParticles; j++) {
Vec3 pos = atomCoordinates[site.getParticle(j)];
origin += pos*originWeights[j];
xdir += pos*xWeights[j];
ydir += pos*yWeights[j];
}
Vec3 localPosition = site.getLocalPosition();
Vec3 origin = atomCoordinates[p1]*originWeights[0] + atomCoordinates[p2]*originWeights[1] + atomCoordinates[p3]*originWeights[2];
Vec3 xdir = atomCoordinates[p1]*xWeights[0] + atomCoordinates[p2]*xWeights[1] + atomCoordinates[p3]*xWeights[2];
Vec3 ydir = atomCoordinates[p1]*yWeights[0] + atomCoordinates[p2]*yWeights[1] + atomCoordinates[p3]*yWeights[2];
Vec3 zdir = xdir.cross(ydir);
xdir /= sqrt(xdir.dot(xdir));
zdir /= sqrt(zdir.dot(zdir));
double normXdir = sqrt(xdir.dot(xdir));
double normZdir = sqrt(zdir.dot(zdir));
if (normXdir > 0.0)
xdir /= normXdir;
if (normZdir > 0.0)
zdir /= normZdir;
ydir = zdir.cross(xdir);
atomCoordinates[i] = origin + xdir*localPosition[0] + ydir*localPosition[1] + zdir*localPosition[2];
}
}
}
void ReferenceVirtualSites::distributeForces(const OpenMM::System& system, const vector<OpenMM::Vec3>& atomCoordinates, vector<OpenMM::Vec3>& forces) {
......@@ -133,71 +141,53 @@ void ReferenceVirtualSites::distributeForces(const OpenMM::System& system, const
// A local coordinates site.
const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
Vec3 originWeights = site.getOriginWeights();
Vec3 wx = site.getXWeights();
Vec3 wy = site.getYWeights();
int numParticles = site.getNumParticles();
vector<double> originWeights, wx, wy;
site.getOriginWeights(originWeights);
site.getXWeights(wx);
site.getYWeights(wy);
Vec3 xdir, ydir;
for (int j = 0; j < numParticles; j++) {
Vec3 pos = atomCoordinates[site.getParticle(j)];
xdir += pos*wx[j];
ydir += pos*wy[j];
}
Vec3 localPosition = site.getLocalPosition();
Vec3 xdir = atomCoordinates[p1]*wx[0] + atomCoordinates[p2]*wx[1] + atomCoordinates[p3]*wx[2];
Vec3 ydir = atomCoordinates[p1]*wy[0] + atomCoordinates[p2]*wy[1] + atomCoordinates[p3]*wy[2];
Vec3 zdir = xdir.cross(ydir);
double invNormXdir = 1.0/sqrt(xdir.dot(xdir));
double invNormZdir = 1.0/sqrt(zdir.dot(zdir));
double normXdir = sqrt(xdir.dot(xdir));
double normZdir = sqrt(zdir.dot(zdir));
double invNormXdir = (normXdir > 0.0 ? 1.0/normXdir : 0.0);
double invNormZdir = (normZdir > 0.0 ? 1.0/normZdir : 0.0);
Vec3 dx = xdir*invNormXdir;
Vec3 dz = zdir*invNormZdir;
Vec3 dy = dz.cross(dx);
// The derivatives for this case are very complicated. They were computed with SymPy then simplified by hand.
double t11 = (wx[0]*ydir[0]-wy[0]*xdir[0])*invNormZdir;
double t12 = (wx[0]*ydir[1]-wy[0]*xdir[1])*invNormZdir;
double t13 = (wx[0]*ydir[2]-wy[0]*xdir[2])*invNormZdir;
double t21 = (wx[1]*ydir[0]-wy[1]*xdir[0])*invNormZdir;
double t22 = (wx[1]*ydir[1]-wy[1]*xdir[1])*invNormZdir;
double t23 = (wx[1]*ydir[2]-wy[1]*xdir[2])*invNormZdir;
double t31 = (wx[2]*ydir[0]-wy[2]*xdir[0])*invNormZdir;
double t32 = (wx[2]*ydir[1]-wy[2]*xdir[1])*invNormZdir;
double t33 = (wx[2]*ydir[2]-wy[2]*xdir[2])*invNormZdir;
double sx1 = t13*dz[1]-t12*dz[2];
double sy1 = t11*dz[2]-t13*dz[0];
double sz1 = t12*dz[0]-t11*dz[1];
double sx2 = t23*dz[1]-t22*dz[2];
double sy2 = t21*dz[2]-t23*dz[0];
double sz2 = t22*dz[0]-t21*dz[1];
double sx3 = t33*dz[1]-t32*dz[2];
double sy3 = t31*dz[2]-t33*dz[0];
double sz3 = t32*dz[0]-t31*dz[1];
Vec3 wxScaled = wx*invNormXdir;
vector<double> wxScaled(numParticles);
for (int j = 0; j < numParticles; j++)
wxScaled[j] = wx[j]*invNormXdir;
Vec3 fp1 = localPosition*f[0];
Vec3 fp2 = localPosition*f[1];
Vec3 fp3 = localPosition*f[2];
forces[p1][0] += fp1[0]*wxScaled[0]*(1-dx[0]*dx[0]) + fp1[2]*(dz[0]*sx1 ) + fp1[1]*((-dx[0]*dy[0] )*wxScaled[0] + dy[0]*sx1 - dx[1]*t12 - dx[2]*t13) + f[0]*originWeights[0];
forces[p1][1] += fp1[0]*wxScaled[0]*( -dx[0]*dx[1]) + fp1[2]*(dz[0]*sy1+t13) + fp1[1]*((-dx[1]*dy[0]-dz[2])*wxScaled[0] + dy[0]*sy1 + dx[1]*t11);
forces[p1][2] += fp1[0]*wxScaled[0]*( -dx[0]*dx[2]) + fp1[2]*(dz[0]*sz1-t12) + fp1[1]*((-dx[2]*dy[0]+dz[1])*wxScaled[0] + dy[0]*sz1 + dx[2]*t11);
forces[p2][0] += fp1[0]*wxScaled[1]*(1-dx[0]*dx[0]) + fp1[2]*(dz[0]*sx2 ) + fp1[1]*((-dx[0]*dy[0] )*wxScaled[1] + dy[0]*sx2 - dx[1]*t22 - dx[2]*t23) + f[0]*originWeights[1];
forces[p2][1] += fp1[0]*wxScaled[1]*( -dx[0]*dx[1]) + fp1[2]*(dz[0]*sy2+t23) + fp1[1]*((-dx[1]*dy[0]-dz[2])*wxScaled[1] + dy[0]*sy2 + dx[1]*t21);
forces[p2][2] += fp1[0]*wxScaled[1]*( -dx[0]*dx[2]) + fp1[2]*(dz[0]*sz2-t22) + fp1[1]*((-dx[2]*dy[0]+dz[1])*wxScaled[1] + dy[0]*sz2 + dx[2]*t21);
forces[p3][0] += fp1[0]*wxScaled[2]*(1-dx[0]*dx[0]) + fp1[2]*(dz[0]*sx3 ) + fp1[1]*((-dx[0]*dy[0] )*wxScaled[2] + dy[0]*sx3 - dx[1]*t32 - dx[2]*t33) + f[0]*originWeights[2];
forces[p3][1] += fp1[0]*wxScaled[2]*( -dx[0]*dx[1]) + fp1[2]*(dz[0]*sy3+t33) + fp1[1]*((-dx[1]*dy[0]-dz[2])*wxScaled[2] + dy[0]*sy3 + dx[1]*t31);
forces[p3][2] += fp1[0]*wxScaled[2]*( -dx[0]*dx[2]) + fp1[2]*(dz[0]*sz3-t32) + fp1[1]*((-dx[2]*dy[0]+dz[1])*wxScaled[2] + dy[0]*sz3 + dx[2]*t31);
forces[p1][0] += fp2[0]*wxScaled[0]*( -dx[1]*dx[0]) + fp2[2]*(dz[1]*sx1-t13) - fp2[1]*(( dx[0]*dy[1]-dz[2])*wxScaled[0] - dy[1]*sx1 - dx[0]*t12);
forces[p1][1] += fp2[0]*wxScaled[0]*(1-dx[1]*dx[1]) + fp2[2]*(dz[1]*sy1 ) - fp2[1]*(( dx[1]*dy[1] )*wxScaled[0] - dy[1]*sy1 + dx[0]*t11 + dx[2]*t13) + f[1]*originWeights[0];
forces[p1][2] += fp2[0]*wxScaled[0]*( -dx[1]*dx[2]) + fp2[2]*(dz[1]*sz1+t11) - fp2[1]*(( dx[2]*dy[1]+dz[0])*wxScaled[0] - dy[1]*sz1 - dx[2]*t12);
forces[p2][0] += fp2[0]*wxScaled[1]*( -dx[1]*dx[0]) + fp2[2]*(dz[1]*sx2-t23) - fp2[1]*(( dx[0]*dy[1]-dz[2])*wxScaled[1] - dy[1]*sx2 - dx[0]*t22);
forces[p2][1] += fp2[0]*wxScaled[1]*(1-dx[1]*dx[1]) + fp2[2]*(dz[1]*sy2 ) - fp2[1]*(( dx[1]*dy[1] )*wxScaled[1] - dy[1]*sy2 + dx[0]*t21 + dx[2]*t23) + f[1]*originWeights[1];
forces[p2][2] += fp2[0]*wxScaled[1]*( -dx[1]*dx[2]) + fp2[2]*(dz[1]*sz2+t21) - fp2[1]*(( dx[2]*dy[1]+dz[0])*wxScaled[1] - dy[1]*sz2 - dx[2]*t22);
forces[p3][0] += fp2[0]*wxScaled[2]*( -dx[1]*dx[0]) + fp2[2]*(dz[1]*sx3-t33) - fp2[1]*(( dx[0]*dy[1]-dz[2])*wxScaled[2] - dy[1]*sx3 - dx[0]*t32);
forces[p3][1] += fp2[0]*wxScaled[2]*(1-dx[1]*dx[1]) + fp2[2]*(dz[1]*sy3 ) - fp2[1]*(( dx[1]*dy[1] )*wxScaled[2] - dy[1]*sy3 + dx[0]*t31 + dx[2]*t33) + f[1]*originWeights[2];
forces[p3][2] += fp2[0]*wxScaled[2]*( -dx[1]*dx[2]) + fp2[2]*(dz[1]*sz3+t31) - fp2[1]*(( dx[2]*dy[1]+dz[0])*wxScaled[2] - dy[1]*sz3 - dx[2]*t32);
forces[p1][0] += fp3[0]*wxScaled[0]*( -dx[2]*dx[0]) + fp3[2]*(dz[2]*sx1+t12) + fp3[1]*((-dx[0]*dy[2]-dz[1])*wxScaled[0] + dy[2]*sx1 + dx[0]*t13);
forces[p1][1] += fp3[0]*wxScaled[0]*( -dx[2]*dx[1]) + fp3[2]*(dz[2]*sy1-t11) + fp3[1]*((-dx[1]*dy[2]+dz[0])*wxScaled[0] + dy[2]*sy1 + dx[1]*t13);
forces[p1][2] += fp3[0]*wxScaled[0]*(1-dx[2]*dx[2]) + fp3[2]*(dz[2]*sz1 ) + fp3[1]*((-dx[2]*dy[2] )*wxScaled[0] + dy[2]*sz1 - dx[0]*t11 - dx[1]*t12) + f[2]*originWeights[0];
forces[p2][0] += fp3[0]*wxScaled[1]*( -dx[2]*dx[0]) + fp3[2]*(dz[2]*sx2+t22) + fp3[1]*((-dx[0]*dy[2]-dz[1])*wxScaled[1] + dy[2]*sx2 + dx[0]*t23);
forces[p2][1] += fp3[0]*wxScaled[1]*( -dx[2]*dx[1]) + fp3[2]*(dz[2]*sy2-t21) + fp3[1]*((-dx[1]*dy[2]+dz[0])*wxScaled[1] + dy[2]*sy2 + dx[1]*t23);
forces[p2][2] += fp3[0]*wxScaled[1]*(1-dx[2]*dx[2]) + fp3[2]*(dz[2]*sz2 ) + fp3[1]*((-dx[2]*dy[2] )*wxScaled[1] + dy[2]*sz2 - dx[0]*t21 - dx[1]*t22) + f[2]*originWeights[1];
forces[p3][0] += fp3[0]*wxScaled[2]*( -dx[2]*dx[0]) + fp3[2]*(dz[2]*sx3+t32) + fp3[1]*((-dx[0]*dy[2]-dz[1])*wxScaled[2] + dy[2]*sx3 + dx[0]*t33);
forces[p3][1] += fp3[0]*wxScaled[2]*( -dx[2]*dx[1]) + fp3[2]*(dz[2]*sy3-t31) + fp3[1]*((-dx[1]*dy[2]+dz[0])*wxScaled[2] + dy[2]*sy3 + dx[1]*t33);
forces[p3][2] += fp3[0]*wxScaled[2]*(1-dx[2]*dx[2]) + fp3[2]*(dz[2]*sz3 ) + fp3[1]*((-dx[2]*dy[2] )*wxScaled[2] + dy[2]*sz3 - dx[0]*t31 - dx[1]*t32) + f[2]*originWeights[2];
for (int j = 0; j < numParticles; j++) {
double t1 = (wx[j]*ydir[0]-wy[j]*xdir[0])*invNormZdir;
double t2 = (wx[j]*ydir[1]-wy[j]*xdir[1])*invNormZdir;
double t3 = (wx[j]*ydir[2]-wy[j]*xdir[2])*invNormZdir;
double sx = t3*dz[1]-t2*dz[2];
double sy = t1*dz[2]-t3*dz[0];
double sz = t2*dz[0]-t1*dz[1];
int p = site.getParticle(j);
forces[p][0] += fp1[0]*wxScaled[j]*(1-dx[0]*dx[0]) + fp1[2]*(dz[0]*sx ) + fp1[1]*((-dx[0]*dy[0] )*wxScaled[j] + dy[0]*sx - dx[1]*t2 - dx[2]*t3) + f[0]*originWeights[j];
forces[p][1] += fp1[0]*wxScaled[j]*( -dx[0]*dx[1]) + fp1[2]*(dz[0]*sy+t3) + fp1[1]*((-dx[1]*dy[0]-dz[2])*wxScaled[j] + dy[0]*sy + dx[1]*t1);
forces[p][2] += fp1[0]*wxScaled[j]*( -dx[0]*dx[2]) + fp1[2]*(dz[0]*sz-t2) + fp1[1]*((-dx[2]*dy[0]+dz[1])*wxScaled[j] + dy[0]*sz + dx[2]*t1);
forces[p][0] += fp2[0]*wxScaled[j]*( -dx[1]*dx[0]) + fp2[2]*(dz[1]*sx-t3) - fp2[1]*(( dx[0]*dy[1]-dz[2])*wxScaled[j] - dy[1]*sx - dx[0]*t2);
forces[p][1] += fp2[0]*wxScaled[j]*(1-dx[1]*dx[1]) + fp2[2]*(dz[1]*sy ) - fp2[1]*(( dx[1]*dy[1] )*wxScaled[j] - dy[1]*sy + dx[0]*t1 + dx[2]*t3) + f[1]*originWeights[j];
forces[p][2] += fp2[0]*wxScaled[j]*( -dx[1]*dx[2]) + fp2[2]*(dz[1]*sz+t1) - fp2[1]*(( dx[2]*dy[1]+dz[0])*wxScaled[j] - dy[1]*sz - dx[2]*t2);
forces[p][0] += fp3[0]*wxScaled[j]*( -dx[2]*dx[0]) + fp3[2]*(dz[2]*sx+t2) + fp3[1]*((-dx[0]*dy[2]-dz[1])*wxScaled[j] + dy[2]*sx + dx[0]*t3);
forces[p][1] += fp3[0]*wxScaled[j]*( -dx[2]*dx[1]) + fp3[2]*(dz[2]*sy-t1) + fp3[1]*((-dx[1]*dy[2]+dz[0])*wxScaled[j] + dy[2]*sy + dx[1]*t3);
forces[p][2] += fp3[0]*wxScaled[j]*(1-dx[2]*dx[2]) + fp3[2]*(dz[2]*sz ) + fp3[1]*((-dx[2]*dy[2] )*wxScaled[j] + dy[2]*sz - dx[0]*t1 - dx[1]*t2) + f[2]*originWeights[j];
}
}
}
}
serialization/src/SystemProxy.cpp
View file @ 0827ad50
......@@ -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-2015 Stanford University and the Authors. *
* Portions copyright (c) 2010-2017 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -71,15 +71,23 @@ void SystemProxy::serialize(const void* object, SerializationNode& node) const {
}
else if (typeid(system.getVirtualSite(i)) == typeid(LocalCoordinatesSite)) {
const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
Vec3 wo = site.getOriginWeights();
Vec3 wx = site.getXWeights();
Vec3 wy = site.getYWeights();
int numParticles = site.getNumParticles();
vector<double> wo, wx, wy;
site.getOriginWeights(wo);
site.getXWeights(wx);
site.getYWeights(wy);
Vec3 p = site.getLocalPosition();
particle.createChildNode("LocalCoordinatesSite").setIntProperty("p1", site.getParticle(0)).setIntProperty("p2", site.getParticle(1)).setIntProperty("p3", site.getParticle(2)).
setDoubleProperty("wo1", wo[0]).setDoubleProperty("wo2", wo[1]).setDoubleProperty("wo3", wo[2]).
setDoubleProperty("wx1", wx[0]).setDoubleProperty("wx2", wx[1]).setDoubleProperty("wx3", wx[2]).
setDoubleProperty("wy1", wy[0]).setDoubleProperty("wy2", wy[1]).setDoubleProperty("wy3", wy[2]).
setDoubleProperty("pos1", p[0]).setDoubleProperty("pos2", p[1]).setDoubleProperty("pos3", p[2]);
SerializationNode& siteNode = particle.createChildNode("LocalCoordinatesSite");
siteNode.setDoubleProperty("pos1", p[0]).setDoubleProperty("pos2", p[1]).setDoubleProperty("pos3", p[2]);
for (int j = 0; j < numParticles; j++) {
stringstream ss;
ss << (j+1);
string index = ss.str();
siteNode.setIntProperty("p"+index, site.getParticle(j));
siteNode.setDoubleProperty("wo"+index, wo[j]);
siteNode.setDoubleProperty("wx"+index, wx[j]);
siteNode.setDoubleProperty("wy"+index, wy[j]);
}
}
}
}
......@@ -120,11 +128,21 @@ void* SystemProxy::deserialize(const SerializationNode& node) const {
else if (vsite.getName() == "OutOfPlaneSite")
system->setVirtualSite(i, new OutOfPlaneSite(vsite.getIntProperty("p1"), vsite.getIntProperty("p2"), vsite.getIntProperty("p3"), vsite.getDoubleProperty("w12"), vsite.getDoubleProperty("w13"), vsite.getDoubleProperty("wc")));
else if (vsite.getName() == "LocalCoordinatesSite") {
Vec3 wo(vsite.getDoubleProperty("wo1"), vsite.getDoubleProperty("wo2"), vsite.getDoubleProperty("wo3"));
Vec3 wx(vsite.getDoubleProperty("wx1"), vsite.getDoubleProperty("wx2"), vsite.getDoubleProperty("wx3"));
Vec3 wy(vsite.getDoubleProperty("wy1"), vsite.getDoubleProperty("wy2"), vsite.getDoubleProperty("wy3"));
vector<int> particles;
vector<double> wo, wx, wy;
for (int j = 0; ; j++) {
stringstream ss;
ss << (j+1);
string index = ss.str();
if (!vsite.hasProperty("p"+index))
break;
particles.push_back(vsite.getIntProperty("p"+index));
wo.push_back(vsite.getDoubleProperty("wo"+index));
wx.push_back(vsite.getDoubleProperty("wx"+index));
wy.push_back(vsite.getDoubleProperty("wy"+index));
}
Vec3 p(vsite.getDoubleProperty("pos1"), vsite.getDoubleProperty("pos2"), vsite.getDoubleProperty("pos3"));
system->setVirtualSite(i, new LocalCoordinatesSite(vsite.getIntProperty("p1"), vsite.getIntProperty("p2"), vsite.getIntProperty("p3"), wo, wx, wy, p));
system->setVirtualSite(i, new LocalCoordinatesSite(particles, wo, wx, wy, p));
}
}
}
......
serialization/tests/TestSerializeSystem.cpp
View file @ 0827ad50
......@@ -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-2017 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -82,6 +82,23 @@ void compareSystems(System& system, System& system2) {
ASSERT_EQUAL(0.1, site7.getWeight12());
ASSERT_EQUAL(0.2, site7.getWeight13());
ASSERT_EQUAL(0.5, site7.getWeightCross());
const LocalCoordinatesSite& site8 = dynamic_cast<const LocalCoordinatesSite&>(system2.getVirtualSite(8));
ASSERT_EQUAL(4, site8.getNumParticles());
ASSERT_EQUAL(4, site8.getParticle(0));
ASSERT_EQUAL(3, site8.getParticle(1));
ASSERT_EQUAL(2, site8.getParticle(2));
ASSERT_EQUAL(1, site8.getParticle(3));
ASSERT_EQUAL(Vec3(-0.5, 1.0, 1.5), site8.getLocalPosition());
vector<double> wo, wx, wy;
site8.getOriginWeights(wo);
site8.getXWeights(wx);
site8.getYWeights(wy);
vector<double> woExpected = {0.1, 0.2, 0.3, 0.4};
vector<double> wxExpected = {-1.0, 0.4, 0.4, 0.2};
vector<double> wyExpected = {0.3, 0.7, 0.0, -1.0};
ASSERT_EQUAL_CONTAINERS(woExpected, wo);
ASSERT_EQUAL_CONTAINERS(wxExpected, wx);
ASSERT_EQUAL_CONTAINERS(wyExpected, wy);
ASSERT_EQUAL(system.getNumForces(), system2.getNumForces());
for (int i = 0; i < system.getNumForces(); i++)
ASSERT(typeid(system.getForce(i)) == typeid(system2.getForce(i)))
......@@ -93,7 +110,7 @@ void testSerialization() {
System system;
for (int i = 0; i < 5; i++)
system.addParticle(0.1*i+1);
for (int i = 0; i < 4; i++)
for (int i = 0; i < 5; i++)
system.addParticle(0.0);
system.addConstraint(0, 1, 3.0);
system.addConstraint(1, 2, 2.5);
......@@ -102,6 +119,7 @@ void testSerialization() {
system.setVirtualSite(5, new TwoParticleAverageSite(0, 1, 0.3, 0.7));
system.setVirtualSite(6, new ThreeParticleAverageSite(2, 4, 3, 0.5, 0.2, 0.3));
system.setVirtualSite(7, new OutOfPlaneSite(0, 3, 1, 0.1, 0.2, 0.5));
system.setVirtualSite(8, new LocalCoordinatesSite({4, 3, 2, 1}, {0.1, 0.2, 0.3, 0.4}, {-1.0, 0.4, 0.4, 0.2}, {0.3, 0.7, 0.0, -1.0}, Vec3(-0.5, 1.0, 1.5)));
system.addForce(new HarmonicBondForce());
// Serialize and then deserialize it, then make sure the systems are identical.
......
tests/TestVirtualSites.h
View file @ 0827ad50
......@@ -206,30 +206,51 @@ void testOutOfPlane() {
}
}
Vec3 computeWeightedPosition(const vector<Vec3>& positions, const vector<double>& weights) {
Vec3 sum;
for (int i = 0; i < weights.size(); i++)
sum += positions[i]*weights[i];
return sum;
}
/**
* Test a LocalCoordinatesSite virtual site.
*/
void testLocalCoordinates() {
const Vec3 originWeights(0.2, 0.3, 0.5);
const Vec3 xWeights(-1.0, 0.5, 0.5);
const Vec3 yWeights(0.0, -1.0, 1.0);
void testLocalCoordinates(int numSiteParticles) {
vector<int> particles;
vector<double> originWeights, xWeights, yWeights;
if (numSiteParticles == 2) {
particles = {0, 1};
originWeights = {0.4, 0.6};
xWeights = {-1.0, 1.0};
yWeights = {1.0, -1.0};
}
else if (numSiteParticles == 3) {
particles = {0, 1, 2};
originWeights = {0.2, 0.3, 0.5};
xWeights = {-1.0, 0.5, 0.5};
yWeights = {0.0, -1.0, 1.0};
}
else if (numSiteParticles == 4) {
particles = {0, 1, 2, 3};
originWeights = {0.2, 0.3, 0.1, 0.4};
xWeights = {-1.0, 0.3, 0.3, 0.4};
yWeights = {0.5, 0.5, -0.5, -0.5};
}
const Vec3 localPosition(0.4, 0.3, 0.2);
System system;
system.addParticle(1.0);
system.addParticle(1.0);
for (int i = 0; i < numSiteParticles; i++)
system.addParticle(1.0);
system.addParticle(0.0);
system.setVirtualSite(3, new LocalCoordinatesSite(0, 1, 2, originWeights, xWeights, yWeights, localPosition));
system.setVirtualSite(numSiteParticles, new LocalCoordinatesSite(particles, originWeights, xWeights, yWeights, localPosition));
CustomExternalForce* forceField = new CustomExternalForce("2*x^2+3*y^2+4*z^2");
system.addForce(forceField);
vector<double> params;
for (int i = 0; i < numSiteParticles+1; i++)
forceField->addParticle(0, params);
forceField->addParticle(1, params);
forceField->addParticle(2, params);
forceField->addParticle(3, params);
LangevinIntegrator integrator(300.0, 0.1, 0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4), positions2(4), positions3(4);
vector<Vec3> positions(numSiteParticles+1), positions2(numSiteParticles+1), positions3(numSiteParticles+1);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
for (int i = 0; i < 100; i++) {
......@@ -237,12 +258,12 @@ void testLocalCoordinates() {
Vec3 xdir, ydir, zdir;
do {
for (int j = 0; j < 3; j++)
for (int j = 0; j < numSiteParticles; j++)
positions[j] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
xdir = positions[0]*xWeights[0] + positions[1]*xWeights[1] + positions[2]*xWeights[2];
ydir = positions[0]*yWeights[0] + positions[1]*yWeights[1] + positions[2]*yWeights[2];
xdir = computeWeightedPosition(positions, xWeights);
ydir = computeWeightedPosition(positions, yWeights);;
zdir = xdir.cross(ydir);
if (sqrt(xdir.dot(xdir)) > 0.1 && sqrt(ydir.dot(ydir)) > 0.1 && sqrt(zdir.dot(zdir)) > 0.1)
if (sqrt(xdir.dot(xdir)) > 0.1 && (numSiteParticles == 2 || (sqrt(ydir.dot(ydir)) > 0.1 && sqrt(zdir.dot(zdir)) > 0.1)))
break; // These positions give a reasonable coordinate system.
} while (true);
context.setPositions(positions);
......@@ -252,23 +273,23 @@ void testLocalCoordinates() {
State state = context.getState(State::Positions | State::Forces);
const vector<Vec3>& pos = state.getPositions();
Vec3 origin = pos[0]*originWeights[0] + pos[1]*originWeights[1] + pos[2]*originWeights[2];
Vec3 origin = computeWeightedPosition(pos, originWeights);;
xdir /= sqrt(xdir.dot(xdir));
zdir /= sqrt(zdir.dot(zdir));
ydir = zdir.cross(xdir);
ASSERT_EQUAL_VEC(origin+xdir*localPosition[0]+ydir*localPosition[1]+zdir*localPosition[2], pos[3], 1e-5);
ASSERT_EQUAL_VEC(origin+xdir*localPosition[0]+ydir*localPosition[1]+zdir*localPosition[2], pos[numSiteParticles], 1e-5);
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
double norm = 0.0;
for (int i = 0; i < 3; ++i) {
for (int i = 0; i < numSiteParticles; ++i) {
Vec3 f = state.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-2;
double step = 0.5*delta/norm;
for (int i = 0; i < 3; ++i) {
for (int i = 0; i < numSiteParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = state.getForces()[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
......@@ -469,7 +490,9 @@ int main(int argc, char* argv[]) {
testTwoParticleAverage();
testThreeParticleAverage();
testOutOfPlane();
testLocalCoordinates();
testLocalCoordinates(2);
testLocalCoordinates(3);
testLocalCoordinates(4);
testConservationLaws();
testOverlappingSites();
runPlatformTests();
......
wrappers/generateWrappers.py
View file @ 0827ad50
......@@ -68,7 +68,15 @@ class WrapperGenerator:
def __init__(self, inputDirname, output):
self.skipClasses = ['OpenMM::Vec3', 'OpenMM::XmlSerializer', 'OpenMM::Kernel', 'OpenMM::KernelImpl', 'OpenMM::KernelFactory', 'OpenMM::ContextImpl', 'OpenMM::SerializationNode', 'OpenMM::SerializationProxy']
self.skipMethods = ['OpenMM::Context::getState', 'OpenMM::Platform::loadPluginsFromDirectory', 'OpenMM::Platform::getPluginLoadFailures', 'OpenMM::Context::createCheckpoint', 'OpenMM::Context::loadCheckpoint', 'OpenMM::Context::getMolecules']
self.skipMethods = ['State OpenMM::Context::getState',
'void OpenMM::Context::createCheckpoint',
'void OpenMM::Context::loadCheckpoint',
'const std::vector<std::vector<int> >& OpenMM::Context::getMolecules',
'static std::vector<std::string> OpenMM::Platform::getPluginLoadFailures',
'static std::vector<std::string> OpenMM::Platform::loadPluginsFromDirectory',
'Vec3 OpenMM::LocalCoordinatesSite::getOriginWeights',
'Vec3 OpenMM::LocalCoordinatesSite::getXWeights',
'Vec3 OpenMM::LocalCoordinatesSite::getYWeights']
self.hideClasses = ['Kernel', 'KernelImpl', 'KernelFactory', 'ContextImpl', 'SerializationNode', 'SerializationProxy']
self.nodeByID={}
......@@ -119,9 +127,7 @@ class WrapperGenerator:
for section in findNodes(classNode, "sectiondef", kind="public-static-func")+findNodes(classNode, "sectiondef", kind="public-func"):
for memberNode in findNodes(section, "memberdef", kind="function", prot="public"):
methodDefinition = getText("definition", memberNode)
shortMethodDefinition = stripOpenMMPrefix(methodDefinition)
methodName = shortMethodDefinition.split()[-1]
if className+'::'+methodName in self.skipMethods:
if methodDefinition in self.skipMethods:
continue
methodList.append(memberNode)
return methodList
......
wrappers/python/simtk/openmm/app/forcefield.py
View file @ 0827ad50
......@@ -647,10 +647,15 @@ class ForceField(object):
numAtoms = 3
self.weights = [float(attrib['weight12']), float(attrib['weight13']), float(attrib['weightCross'])]
elif self.type == 'localCoords':
numAtoms = 3
self.originWeights = [float(attrib['wo1']), float(attrib['wo2']), float(attrib['wo3'])]
self.xWeights = [float(attrib['wx1']), float(attrib['wx2']), float(attrib['wx3'])]
self.yWeights = [float(attrib['wy1']), float(attrib['wy2']), float(attrib['wy3'])]
numAtoms = 0
self.originWeights = []
self.xWeights = []
self.yWeights = []
while ('wo%d' % (numAtoms+1)) in attrib:
numAtoms += 1
self.originWeights.append(float(attrib['wo%d' % numAtoms]))
self.xWeights.append(float(attrib['wx%d' % numAtoms]))
self.yWeights.append(float(attrib['wy%d' % numAtoms]))
self.localPos = [float(attrib['p1']), float(attrib['p2']), float(attrib['p3'])]
else:
raise ValueError('Unknown virtual site type: %s' % self.type)
......@@ -1255,11 +1260,7 @@ class ForceField(object):
elif site.type == 'outOfPlane':
sys.setVirtualSite(index, mm.OutOfPlaneSite(atoms[0], atoms[1], atoms[2], site.weights[0], site.weights[1], site.weights[2]))
elif site.type == 'localCoords':
sys.setVirtualSite(index, mm.LocalCoordinatesSite(atoms[0], atoms[1], atoms[2],
mm.Vec3(site.originWeights[0], site.originWeights[1], site.originWeights[2]),
mm.Vec3(site.xWeights[0], site.xWeights[1], site.xWeights[2]),
mm.Vec3(site.yWeights[0], site.yWeights[1], site.yWeights[2]),
mm.Vec3(site.localPos[0], site.localPos[1], site.localPos[2])))
sys.setVirtualSite(index, mm.LocalCoordinatesSite(atoms, site.originWeights, site.xWeights, site.yWeights, site.localPos))
# Add forces to the System
......
wrappers/python/src/swig_doxygen/swigInputConfig.py
View file @ 0827ad50
......@@ -108,6 +108,9 @@ SKIP_METHODS = [('State', 'getPositions'),
('IntegrateDrudeSCFStepKernel',),
('XmlSerializer', 'serialize'),
('XmlSerializer', 'deserialize'),
('LocalCoordinatesSite', 'getOriginWeights', 0),
('LocalCoordinatesSite', 'getXWeights', 0),
('LocalCoordinatesSite', 'getYWeights', 0),
]
# The build script assumes method args that are non-const references are
......
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