Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
tsoc
openmm
Commits
08e8b206
Commit
08e8b206
authored
Mar 31, 2017
by
peastman
Committed by
GitHub
Mar 31, 2017
Browse files
Merge pull request #1769 from peastman/loops
Use C++11 style loops
parents
55ee0b9f
083bc501
Changes
96
Expand all
Hide whitespace changes
Inline
Side-by-side
Showing
20 changed files
with
237 additions
and
293 deletions
+237
-293
platforms/reference/src/SimTKReference/ReferenceCustomAngleIxn.cpp
.../reference/src/SimTKReference/ReferenceCustomAngleIxn.cpp
+4
-4
platforms/reference/src/SimTKReference/ReferenceCustomBondIxn.cpp
...s/reference/src/SimTKReference/ReferenceCustomBondIxn.cpp
+4
-4
platforms/reference/src/SimTKReference/ReferenceCustomCentroidBondIxn.cpp
...nce/src/SimTKReference/ReferenceCustomCentroidBondIxn.cpp
+16
-26
platforms/reference/src/SimTKReference/ReferenceCustomCompoundBondIxn.cpp
...nce/src/SimTKReference/ReferenceCustomCompoundBondIxn.cpp
+16
-26
platforms/reference/src/SimTKReference/ReferenceCustomDynamics.cpp
.../reference/src/SimTKReference/ReferenceCustomDynamics.cpp
+11
-11
platforms/reference/src/SimTKReference/ReferenceCustomExternalIxn.cpp
...ference/src/SimTKReference/ReferenceCustomExternalIxn.cpp
+10
-10
platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.cpp
...rms/reference/src/SimTKReference/ReferenceCustomGBIxn.cpp
+11
-14
platforms/reference/src/SimTKReference/ReferenceCustomHbondIxn.cpp
.../reference/src/SimTKReference/ReferenceCustomHbondIxn.cpp
+6
-6
platforms/reference/src/SimTKReference/ReferenceCustomManyParticleIxn.cpp
...nce/src/SimTKReference/ReferenceCustomManyParticleIxn.cpp
+16
-26
platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp
...erence/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp
+9
-10
platforms/reference/src/SimTKReference/ReferenceCustomTorsionIxn.cpp
...eference/src/SimTKReference/ReferenceCustomTorsionIxn.cpp
+4
-4
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
...ms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
+5
-8
platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp
...erence/src/SimTKReference/ReferenceMonteCarloBarostat.cpp
+6
-6
platforms/reference/src/SimTKReference/ReferenceNeighborList.cpp
...ms/reference/src/SimTKReference/ReferenceNeighborList.cpp
+3
-3
plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp
plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp
+7
-7
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
+26
-25
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaExtrapolatedPolarization.cpp
...rms/cuda/tests/TestCudaAmoebaExtrapolatedPolarization.cpp
+2
-2
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp
...rms/cuda/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp
+2
-2
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaMultipoleForce.cpp
...eba/platforms/cuda/tests/TestCudaAmoebaMultipoleForce.cpp
+2
-2
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp
...ence/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp
+77
-97
No files found.
platforms/reference/src/SimTKReference/ReferenceCustomAngleIxn.cpp
View file @
08e8b206
...
...
@@ -47,10 +47,10 @@ ReferenceCustomAngleIxn::ReferenceCustomAngleIxn(const Lepton::CompiledExpressio
expressionSet
.
registerExpression
(
this
->
energyParamDerivExpressions
[
i
]);
thetaIndex
=
expressionSet
.
getVariableIndex
(
"theta"
);
numParameters
=
parameterNames
.
size
();
for
(
int
i
=
0
;
i
<
(
int
)
numP
arameter
s
;
i
++
)
angleParamIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
param
eterNames
[
i
]
));
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globalParameters
.
begin
();
iter
!=
globalParameters
.
end
();
++
iter
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
auto
&
param
:
p
arameter
Names
)
angleParamIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
param
));
for
(
auto
&
param
:
globalParameters
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
param
.
first
),
param
.
second
);
}
/**---------------------------------------------------------------------------------------
...
...
platforms/reference/src/SimTKReference/ReferenceCustomBondIxn.cpp
View file @
08e8b206
...
...
@@ -48,10 +48,10 @@ ReferenceCustomBondIxn::ReferenceCustomBondIxn(const Lepton::CompiledExpression&
expressionSet
.
registerExpression
(
this
->
energyParamDerivExpressions
[
i
]);
rIndex
=
expressionSet
.
getVariableIndex
(
"r"
);
numParameters
=
parameterNames
.
size
();
for
(
int
i
=
0
;
i
<
(
int
)
numP
arameter
s
;
i
++
)
bondParamIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
param
eterNames
[
i
]
));
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globalParameters
.
begin
();
iter
!=
globalParameters
.
end
();
++
iter
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
auto
&
param
:
p
arameter
Names
)
bondParamIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
param
));
for
(
auto
&
param
:
globalParameters
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
param
.
first
),
param
.
second
);
}
/**---------------------------------------------------------------------------------------
...
...
platforms/reference/src/SimTKReference/ReferenceCustomCentroidBondIxn.cpp
View file @
08e8b206
...
...
@@ -56,12 +56,12 @@ ReferenceCustomCentroidBondIxn::ReferenceCustomCentroidBondIxn(int numGroupsPerB
positionTerms
.
push_back
(
ReferenceCustomCentroidBondIxn
::
PositionTermInfo
(
yname
.
str
(),
i
,
1
,
energyExpression
.
differentiate
(
yname
.
str
()).
createCompiledExpression
()));
positionTerms
.
push_back
(
ReferenceCustomCentroidBondIxn
::
PositionTermInfo
(
zname
.
str
(),
i
,
2
,
energyExpression
.
differentiate
(
zname
.
str
()).
createCompiledExpression
()));
}
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
distances
.
begin
();
iter
!=
distances
.
end
();
++
iter
)
distanceTerms
.
push_back
(
ReferenceCustomCentroidBondIxn
::
DistanceTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpression
.
differentiate
(
i
ter
->
first
).
createCompiledExpression
()));
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
angles
.
begin
();
iter
!=
angles
.
end
();
++
iter
)
angleTerms
.
push_back
(
ReferenceCustomCentroidBondIxn
::
AngleTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpression
.
differentiate
(
i
ter
->
first
).
createCompiledExpression
()));
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
dihedrals
.
begin
();
iter
!=
dihedrals
.
end
();
++
iter
)
dihedralTerms
.
push_back
(
ReferenceCustomCentroidBondIxn
::
DihedralTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpression
.
differentiate
(
i
ter
->
first
).
createCompiledExpression
()));
for
(
auto
&
term
:
distances
)
distanceTerms
.
push_back
(
ReferenceCustomCentroidBondIxn
::
DistanceTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpression
.
differentiate
(
ter
m
.
first
).
createCompiledExpression
()));
for
(
auto
&
term
:
angles
)
angleTerms
.
push_back
(
ReferenceCustomCentroidBondIxn
::
AngleTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpression
.
differentiate
(
ter
m
.
first
).
createCompiledExpression
()));
for
(
auto
&
term
:
dihedrals
)
dihedralTerms
.
push_back
(
ReferenceCustomCentroidBondIxn
::
DihedralTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpression
.
differentiate
(
ter
m
.
first
).
createCompiledExpression
()));
for
(
int
i
=
0
;
i
<
positionTerms
.
size
();
i
++
)
{
expressionSet
.
registerExpression
(
positionTerms
[
i
].
forceExpression
);
positionTerms
[
i
].
index
=
expressionSet
.
getVariableIndex
(
positionTerms
[
i
].
name
);
...
...
@@ -108,8 +108,8 @@ void ReferenceCustomCentroidBondIxn::calculatePairIxn(vector<Vec3>& atomCoordina
// Compute the forces on groups.
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globalParameters
.
begin
();
iter
!=
globalParameters
.
end
();
++
iter
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
auto
&
param
:
globalParameters
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
param
.
first
),
param
.
second
);
vector
<
Vec3
>
groupForces
(
numGroups
);
int
numBonds
=
bondGroups
.
size
();
for
(
int
bond
=
0
;
bond
<
numBonds
;
bond
++
)
{
...
...
@@ -131,23 +131,18 @@ void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<Vec3>& gro
// Compute all of the variables the energy can depend on.
const
vector
<
int
>&
groups
=
bondGroups
[
bond
];
for
(
int
i
=
0
;
i
<
(
int
)
positionTerms
.
size
();
i
++
)
{
const
PositionTermInfo
&
term
=
positionTerms
[
i
];
for
(
auto
&
term
:
positionTerms
)
expressionSet
.
setVariable
(
term
.
index
,
groupCenters
[
groups
[
term
.
group
]][
term
.
component
]);
}
for
(
int
i
=
0
;
i
<
(
int
)
distanceTerms
.
size
();
i
++
)
{
const
DistanceTermInfo
&
term
=
distanceTerms
[
i
];
for
(
auto
&
term
:
distanceTerms
)
{
computeDelta
(
groups
[
term
.
g1
],
groups
[
term
.
g2
],
term
.
delta
,
groupCenters
);
expressionSet
.
setVariable
(
term
.
index
,
term
.
delta
[
ReferenceForce
::
RIndex
]);
}
for
(
int
i
=
0
;
i
<
(
int
)
angleTerms
.
size
();
i
++
)
{
const
AngleTermInfo
&
term
=
angleTerms
[
i
];
for
(
auto
&
term
:
angleTerms
)
{
computeDelta
(
groups
[
term
.
g1
],
groups
[
term
.
g2
],
term
.
delta1
,
groupCenters
);
computeDelta
(
groups
[
term
.
g3
],
groups
[
term
.
g2
],
term
.
delta2
,
groupCenters
);
expressionSet
.
setVariable
(
term
.
index
,
computeAngle
(
term
.
delta1
,
term
.
delta2
));
}
for
(
int
i
=
0
;
i
<
(
int
)
dihedralTerms
.
size
();
i
++
)
{
const
DihedralTermInfo
&
term
=
dihedralTerms
[
i
];
for
(
auto
&
term
:
dihedralTerms
)
{
computeDelta
(
groups
[
term
.
g2
],
groups
[
term
.
g1
],
term
.
delta1
,
groupCenters
);
computeDelta
(
groups
[
term
.
g2
],
groups
[
term
.
g3
],
term
.
delta2
,
groupCenters
);
computeDelta
(
groups
[
term
.
g4
],
groups
[
term
.
g3
],
term
.
delta3
,
groupCenters
);
...
...
@@ -158,15 +153,12 @@ void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<Vec3>& gro
// Apply forces based on individual particle coordinates.
for
(
int
i
=
0
;
i
<
(
int
)
positionTerms
.
size
();
i
++
)
{
const
PositionTermInfo
&
term
=
positionTerms
[
i
];
for
(
auto
&
term
:
positionTerms
)
forces
[
groups
[
term
.
group
]][
term
.
component
]
-=
term
.
forceExpression
.
evaluate
();
}
// Apply forces based on distances.
for
(
int
i
=
0
;
i
<
(
int
)
distanceTerms
.
size
();
i
++
)
{
const
DistanceTermInfo
&
term
=
distanceTerms
[
i
];
for
(
auto
&
term
:
distanceTerms
)
{
double
dEdR
=
term
.
forceExpression
.
evaluate
()
/
(
term
.
delta
[
ReferenceForce
::
RIndex
]);
for
(
int
i
=
0
;
i
<
3
;
i
++
)
{
double
force
=
-
dEdR
*
term
.
delta
[
i
];
...
...
@@ -177,8 +169,7 @@ void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<Vec3>& gro
// Apply forces based on angles.
for
(
int
i
=
0
;
i
<
(
int
)
angleTerms
.
size
();
i
++
)
{
const
AngleTermInfo
&
term
=
angleTerms
[
i
];
for
(
auto
&
term
:
angleTerms
)
{
double
dEdTheta
=
term
.
forceExpression
.
evaluate
();
double
thetaCross
[
ReferenceForce
::
LastDeltaRIndex
];
SimTKOpenMMUtilities
::
crossProductVector3
(
term
.
delta1
,
term
.
delta2
,
thetaCross
);
...
...
@@ -204,8 +195,7 @@ void ReferenceCustomCentroidBondIxn::calculateOneIxn(int bond, vector<Vec3>& gro
// Apply forces based on dihedrals.
for
(
int
i
=
0
;
i
<
(
int
)
dihedralTerms
.
size
();
i
++
)
{
const
DihedralTermInfo
&
term
=
dihedralTerms
[
i
];
for
(
auto
&
term
:
dihedralTerms
)
{
double
dEdTheta
=
term
.
forceExpression
.
evaluate
();
double
internalF
[
4
][
3
];
double
forceFactors
[
4
];
...
...
platforms/reference/src/SimTKReference/ReferenceCustomCompoundBondIxn.cpp
View file @
08e8b206
...
...
@@ -61,12 +61,12 @@ ReferenceCustomCompoundBondIxn::ReferenceCustomCompoundBondIxn(int numParticlesP
particleTerms
.
push_back
(
ReferenceCustomCompoundBondIxn
::
ParticleTermInfo
(
yname
.
str
(),
i
,
1
,
energyExpression
.
differentiate
(
yname
.
str
()).
createCompiledExpression
()));
particleTerms
.
push_back
(
ReferenceCustomCompoundBondIxn
::
ParticleTermInfo
(
zname
.
str
(),
i
,
2
,
energyExpression
.
differentiate
(
zname
.
str
()).
createCompiledExpression
()));
}
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
distances
.
begin
();
iter
!=
distances
.
end
();
++
iter
)
distanceTerms
.
push_back
(
ReferenceCustomCompoundBondIxn
::
DistanceTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpression
.
differentiate
(
i
ter
->
first
).
createCompiledExpression
()));
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
angles
.
begin
();
iter
!=
angles
.
end
();
++
iter
)
angleTerms
.
push_back
(
ReferenceCustomCompoundBondIxn
::
AngleTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpression
.
differentiate
(
i
ter
->
first
).
createCompiledExpression
()));
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
dihedrals
.
begin
();
iter
!=
dihedrals
.
end
();
++
iter
)
dihedralTerms
.
push_back
(
ReferenceCustomCompoundBondIxn
::
DihedralTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpression
.
differentiate
(
i
ter
->
first
).
createCompiledExpression
()));
for
(
auto
&
term
:
distances
)
distanceTerms
.
push_back
(
ReferenceCustomCompoundBondIxn
::
DistanceTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpression
.
differentiate
(
ter
m
.
first
).
createCompiledExpression
()));
for
(
auto
&
term
:
angles
)
angleTerms
.
push_back
(
ReferenceCustomCompoundBondIxn
::
AngleTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpression
.
differentiate
(
ter
m
.
first
).
createCompiledExpression
()));
for
(
auto
&
term
:
dihedrals
)
dihedralTerms
.
push_back
(
ReferenceCustomCompoundBondIxn
::
DihedralTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpression
.
differentiate
(
ter
m
.
first
).
createCompiledExpression
()));
for
(
int
i
=
0
;
i
<
particleTerms
.
size
();
i
++
)
{
expressionSet
.
registerExpression
(
particleTerms
[
i
].
forceExpression
);
particleTerms
[
i
].
index
=
expressionSet
.
getVariableIndex
(
particleTerms
[
i
].
name
);
...
...
@@ -119,8 +119,8 @@ void ReferenceCustomCompoundBondIxn::setPeriodic(OpenMM::Vec3* vectors) {
void
ReferenceCustomCompoundBondIxn
::
calculatePairIxn
(
vector
<
Vec3
>&
atomCoordinates
,
double
**
bondParameters
,
const
map
<
string
,
double
>&
globalParameters
,
vector
<
Vec3
>&
forces
,
double
*
totalEnergy
,
double
*
energyParamDerivs
)
{
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globalParameters
.
begin
();
iter
!=
globalParameters
.
end
();
++
iter
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
auto
&
param
:
globalParameters
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
param
.
first
),
param
.
second
);
int
numBonds
=
bondAtoms
.
size
();
for
(
int
bond
=
0
;
bond
<
numBonds
;
bond
++
)
{
for
(
int
i
=
0
;
i
<
numParameters
;
i
++
)
...
...
@@ -146,23 +146,18 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<Vec3>& ato
// Compute all of the variables the energy can depend on.
const
vector
<
int
>&
atoms
=
bondAtoms
[
bond
];
for
(
int
i
=
0
;
i
<
(
int
)
particleTerms
.
size
();
i
++
)
{
const
ParticleTermInfo
&
term
=
particleTerms
[
i
];
for
(
auto
&
term
:
particleTerms
)
expressionSet
.
setVariable
(
term
.
index
,
atomCoordinates
[
atoms
[
term
.
atom
]][
term
.
component
]);
}
for
(
int
i
=
0
;
i
<
(
int
)
distanceTerms
.
size
();
i
++
)
{
const
DistanceTermInfo
&
term
=
distanceTerms
[
i
];
for
(
auto
&
term
:
distanceTerms
)
{
computeDelta
(
atoms
[
term
.
p1
],
atoms
[
term
.
p2
],
term
.
delta
,
atomCoordinates
);
expressionSet
.
setVariable
(
term
.
index
,
term
.
delta
[
ReferenceForce
::
RIndex
]);
}
for
(
int
i
=
0
;
i
<
(
int
)
angleTerms
.
size
();
i
++
)
{
const
AngleTermInfo
&
term
=
angleTerms
[
i
];
for
(
auto
&
term
:
angleTerms
)
{
computeDelta
(
atoms
[
term
.
p1
],
atoms
[
term
.
p2
],
term
.
delta1
,
atomCoordinates
);
computeDelta
(
atoms
[
term
.
p3
],
atoms
[
term
.
p2
],
term
.
delta2
,
atomCoordinates
);
expressionSet
.
setVariable
(
term
.
index
,
computeAngle
(
term
.
delta1
,
term
.
delta2
));
}
for
(
int
i
=
0
;
i
<
(
int
)
dihedralTerms
.
size
();
i
++
)
{
const
DihedralTermInfo
&
term
=
dihedralTerms
[
i
];
for
(
auto
&
term
:
dihedralTerms
)
{
computeDelta
(
atoms
[
term
.
p2
],
atoms
[
term
.
p1
],
term
.
delta1
,
atomCoordinates
);
computeDelta
(
atoms
[
term
.
p2
],
atoms
[
term
.
p3
],
term
.
delta2
,
atomCoordinates
);
computeDelta
(
atoms
[
term
.
p4
],
atoms
[
term
.
p3
],
term
.
delta3
,
atomCoordinates
);
...
...
@@ -173,15 +168,12 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<Vec3>& ato
// Apply forces based on individual particle coordinates.
for
(
int
i
=
0
;
i
<
(
int
)
particleTerms
.
size
();
i
++
)
{
const
ParticleTermInfo
&
term
=
particleTerms
[
i
];
for
(
auto
&
term
:
particleTerms
)
forces
[
atoms
[
term
.
atom
]][
term
.
component
]
-=
term
.
forceExpression
.
evaluate
();
}
// Apply forces based on distances.
for
(
int
i
=
0
;
i
<
(
int
)
distanceTerms
.
size
();
i
++
)
{
const
DistanceTermInfo
&
term
=
distanceTerms
[
i
];
for
(
auto
&
term
:
distanceTerms
)
{
double
dEdR
=
term
.
forceExpression
.
evaluate
()
/
(
term
.
delta
[
ReferenceForce
::
RIndex
]);
for
(
int
i
=
0
;
i
<
3
;
i
++
)
{
double
force
=
-
dEdR
*
term
.
delta
[
i
];
...
...
@@ -192,8 +184,7 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<Vec3>& ato
// Apply forces based on angles.
for
(
int
i
=
0
;
i
<
(
int
)
angleTerms
.
size
();
i
++
)
{
const
AngleTermInfo
&
term
=
angleTerms
[
i
];
for
(
auto
&
term
:
angleTerms
)
{
double
dEdTheta
=
term
.
forceExpression
.
evaluate
();
double
thetaCross
[
ReferenceForce
::
LastDeltaRIndex
];
SimTKOpenMMUtilities
::
crossProductVector3
(
term
.
delta1
,
term
.
delta2
,
thetaCross
);
...
...
@@ -219,8 +210,7 @@ void ReferenceCustomCompoundBondIxn::calculateOneIxn(int bond, vector<Vec3>& ato
// Apply forces based on dihedrals.
for
(
int
i
=
0
;
i
<
(
int
)
dihedralTerms
.
size
();
i
++
)
{
const
DihedralTermInfo
&
term
=
dihedralTerms
[
i
];
for
(
auto
&
term
:
dihedralTerms
)
{
double
dEdTheta
=
term
.
forceExpression
.
evaluate
();
double
internalF
[
4
][
3
];
double
forceFactors
[
4
];
...
...
platforms/reference/src/SimTKReference/ReferenceCustomDynamics.cpp
View file @
08e8b206
...
...
@@ -175,8 +175,8 @@ ExpressionTreeNode ReferenceCustomDynamics::replaceDerivFunctions(const Expressi
}
else
{
vector
<
ExpressionTreeNode
>
children
;
for
(
int
i
=
0
;
i
<
(
int
)
node
.
getChildren
()
.
size
();
i
++
)
children
.
push_back
(
replaceDerivFunctions
(
node
.
getChildren
()[
i
]
,
context
));
for
(
auto
&
child
:
node
.
getChildren
())
children
.
push_back
(
replaceDerivFunctions
(
child
,
context
));
return
ExpressionTreeNode
(
op
.
clone
(),
children
);
}
}
...
...
@@ -204,8 +204,8 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
initialize
(
context
,
masses
,
globals
);
int
numSteps
=
stepType
.
size
();
globals
.
insert
(
context
.
getParameters
().
begin
(),
context
.
getParameters
().
end
());
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globals
.
begin
();
iter
!=
globals
.
end
();
++
iter
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
auto
&
global
:
globals
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
global
.
first
),
global
.
second
);
oldPos
=
atomCoordinates
;
// Loop over steps and execute them.
...
...
@@ -276,8 +276,8 @@ void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, ve
recordChangedParameters
(
context
,
globals
);
context
.
updateContextState
();
globals
.
insert
(
context
.
getParameters
().
begin
(),
context
.
getParameters
().
end
());
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globals
.
begin
();
iter
!=
globals
.
end
();
++
iter
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
auto
&
global
:
globals
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
global
.
first
),
global
.
second
);
break
;
}
case
CustomIntegrator
::
IfBlockStart
:
{
...
...
@@ -355,10 +355,10 @@ bool ReferenceCustomDynamics::evaluateCondition(int step) {
* Check which context parameters have changed and register them with the context.
*/
void
ReferenceCustomDynamics
::
recordChangedParameters
(
OpenMM
::
ContextImpl
&
context
,
std
::
map
<
std
::
string
,
double
>&
globals
)
{
for
(
map
<
string
,
double
>::
const_iterator
iter
=
context
.
getParameters
().
begin
();
iter
!=
context
.
getParameters
()
.
end
();
++
iter
)
{
string
name
=
iter
->
first
;
for
(
auto
&
param
:
context
.
getParameters
())
{
string
name
=
param
.
first
;
double
value
=
globals
[
name
];
if
(
value
!=
iter
->
second
)
if
(
value
!=
param
.
second
)
context
.
setParameter
(
name
,
globals
[
name
]);
}
}
...
...
@@ -385,8 +385,8 @@ double ReferenceCustomDynamics::computeKineticEnergy(OpenMM::ContextImpl& contex
if
(
invalidatesForces
.
size
()
==
0
)
initialize
(
context
,
masses
,
globals
);
globals
.
insert
(
context
.
getParameters
().
begin
(),
context
.
getParameters
().
end
());
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globals
.
begin
();
iter
!=
globals
.
end
();
++
iter
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
auto
&
global
:
globals
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
global
.
first
),
global
.
second
);
if
(
kineticEnergyNeedsForce
)
{
energy
=
context
.
calcForcesAndEnergy
(
true
,
true
,
-
1
);
forcesAreValid
=
true
;
...
...
platforms/reference/src/SimTKReference/ReferenceCustomExternalIxn.cpp
View file @
08e8b206
...
...
@@ -57,17 +57,17 @@ ReferenceCustomExternalIxn::ReferenceCustomExternalIxn(const Lepton::CompiledExp
forceZY
=
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionZ
,
"y"
);
forceZZ
=
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionZ
,
"z"
);
numParameters
=
parameterNames
.
size
();
for
(
int
i
=
0
;
i
<
(
int
)
numP
arameter
s
;
i
++
)
{
energyParams
.
push_back
(
ReferenceForce
::
getVariablePointer
(
this
->
energyExpression
,
param
eterNames
[
i
]
));
forceXParams
.
push_back
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionX
,
param
eterNames
[
i
]
));
forceYParams
.
push_back
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionY
,
param
eterNames
[
i
]
));
forceZParams
.
push_back
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionZ
,
param
eterNames
[
i
]
));
for
(
auto
&
param
:
p
arameter
Names
)
{
energyParams
.
push_back
(
ReferenceForce
::
getVariablePointer
(
this
->
energyExpression
,
param
));
forceXParams
.
push_back
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionX
,
param
));
forceYParams
.
push_back
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionY
,
param
));
forceZParams
.
push_back
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionZ
,
param
));
}
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globalParameters
.
begin
();
iter
!=
globalParameters
.
end
();
++
iter
)
{
ReferenceForce
::
setVariable
(
ReferenceForce
::
getVariablePointer
(
this
->
energyExpression
,
iter
->
first
),
iter
->
second
);
ReferenceForce
::
setVariable
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionX
,
iter
->
first
),
iter
->
second
);
ReferenceForce
::
setVariable
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionY
,
iter
->
first
),
iter
->
second
);
ReferenceForce
::
setVariable
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionZ
,
iter
->
first
),
iter
->
second
);
for
(
auto
&
param
:
globalParameters
)
{
ReferenceForce
::
setVariable
(
ReferenceForce
::
getVariablePointer
(
this
->
energyExpression
,
param
.
first
),
param
.
second
);
ReferenceForce
::
setVariable
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionX
,
param
.
first
),
param
.
second
);
ReferenceForce
::
setVariable
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionY
,
param
.
first
),
param
.
second
);
ReferenceForce
::
setVariable
(
ReferenceForce
::
getVariablePointer
(
this
->
forceExpressionZ
,
param
.
first
),
param
.
second
);
}
}
...
...
platforms/reference/src/SimTKReference/ReferenceCustomGBIxn.cpp
View file @
08e8b206
...
...
@@ -84,19 +84,19 @@ ReferenceCustomGBIxn::ReferenceCustomGBIxn(const vector<Lepton::CompiledExpressi
xIndex
=
expressionSet
.
getVariableIndex
(
"x"
);
yIndex
=
expressionSet
.
getVariableIndex
(
"y"
);
zIndex
=
expressionSet
.
getVariableIndex
(
"z"
);
for
(
int
i
=
0
;
i
<
(
int
)
parameterNames
.
size
();
i
++
)
{
paramIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
param
eterNames
[
i
]
));
for
(
auto
&
param
:
parameterNames
)
{
paramIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
param
));
for
(
int
j
=
1
;
j
<
3
;
j
++
)
{
stringstream
name
;
name
<<
param
eterNames
[
i
]
<<
j
;
name
<<
param
<<
j
;
particleParamIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
name
.
str
()));
}
}
for
(
int
i
=
0
;
i
<
(
int
)
valueNames
.
size
();
i
++
)
{
valueIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
value
Names
[
i
]
));
for
(
auto
&
value
:
valueNames
)
{
valueIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
value
));
for
(
int
j
=
1
;
j
<
3
;
j
++
)
{
stringstream
name
;
name
<<
value
Names
[
i
]
<<
j
;
name
<<
value
<<
j
;
particleValueIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
name
.
str
()));
}
}
...
...
@@ -153,8 +153,8 @@ ReferenceCustomGBIxn::~ReferenceCustomGBIxn() {
void
ReferenceCustomGBIxn
::
calculateIxn
(
int
numberOfAtoms
,
vector
<
Vec3
>&
atomCoordinates
,
double
**
atomParameters
,
const
vector
<
set
<
int
>
>&
exclusions
,
map
<
string
,
double
>&
globalParameters
,
vector
<
Vec3
>&
forces
,
double
*
totalEnergy
,
double
*
energyParamDerivs
)
{
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globalParameters
.
begin
();
iter
!=
globalParameters
.
end
();
++
iter
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
auto
&
param
:
globalParameters
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
param
.
first
),
param
.
second
);
// Initialize arrays for storing values.
...
...
@@ -225,8 +225,7 @@ void ReferenceCustomGBIxn::calculateParticlePairValue(int index, int numAtoms, v
if
(
cutoff
)
{
// Loop over all pairs in the neighbor list.
for
(
int
i
=
0
;
i
<
(
int
)
neighborList
->
size
();
i
++
)
{
OpenMM
::
AtomPair
pair
=
(
*
neighborList
)[
i
];
for
(
auto
&
pair
:
*
neighborList
)
{
if
(
useExclusions
&&
exclusions
[
pair
.
first
].
find
(
pair
.
second
)
!=
exclusions
[
pair
.
first
].
end
())
continue
;
calculateOnePairValue
(
index
,
pair
.
first
,
pair
.
second
,
atomCoordinates
,
atomParameters
);
...
...
@@ -306,8 +305,7 @@ void ReferenceCustomGBIxn::calculateParticlePairEnergyTerm(int index, int numAto
if
(
cutoff
)
{
// Loop over all pairs in the neighbor list.
for
(
int
i
=
0
;
i
<
(
int
)
neighborList
->
size
();
i
++
)
{
OpenMM
::
AtomPair
pair
=
(
*
neighborList
)[
i
];
for
(
auto
&
pair
:
*
neighborList
)
{
if
(
useExclusions
&&
exclusions
[
pair
.
first
].
find
(
pair
.
second
)
!=
exclusions
[
pair
.
first
].
end
())
continue
;
calculateOnePairEnergyTerm
(
index
,
pair
.
first
,
pair
.
second
,
atomCoordinates
,
atomParameters
,
forces
,
totalEnergy
,
energyParamDerivs
);
...
...
@@ -377,8 +375,7 @@ void ReferenceCustomGBIxn::calculateChainRuleForces(int numAtoms, vector<Vec3>&
if
(
cutoff
)
{
// Loop over all pairs in the neighbor list.
for
(
int
i
=
0
;
i
<
(
int
)
neighborList
->
size
();
i
++
)
{
OpenMM
::
AtomPair
pair
=
(
*
neighborList
)[
i
];
for
(
auto
&
pair
:
*
neighborList
)
{
bool
isExcluded
=
(
exclusions
[
pair
.
first
].
find
(
pair
.
second
)
!=
exclusions
[
pair
.
first
].
end
());
calculateOnePairChainRule
(
pair
.
first
,
pair
.
second
,
atomCoordinates
,
atomParameters
,
forces
,
isExcluded
);
calculateOnePairChainRule
(
pair
.
second
,
pair
.
first
,
atomCoordinates
,
atomParameters
,
forces
,
isExcluded
);
...
...
platforms/reference/src/SimTKReference/ReferenceCustomHbondIxn.cpp
View file @
08e8b206
...
...
@@ -49,12 +49,12 @@ ReferenceCustomHbondIxn::ReferenceCustomHbondIxn(const vector<vector<int> >& don
const
map
<
string
,
vector
<
int
>
>&
distances
,
const
map
<
string
,
vector
<
int
>
>&
angles
,
const
map
<
string
,
vector
<
int
>
>&
dihedrals
)
:
cutoff
(
false
),
periodic
(
false
),
donorAtoms
(
donorAtoms
),
acceptorAtoms
(
acceptorAtoms
),
energyExpression
(
energyExpression
.
createProgram
()),
donorParamNames
(
donorParameterNames
),
acceptorParamNames
(
acceptorParameterNames
)
{
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
distances
.
begin
();
iter
!=
distances
.
end
();
++
iter
)
distanceTerms
.
push_back
(
ReferenceCustomHbondIxn
::
DistanceTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpression
.
differentiate
(
i
ter
->
first
).
optimize
().
createProgram
()));
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
angles
.
begin
();
iter
!=
angles
.
end
();
++
iter
)
angleTerms
.
push_back
(
ReferenceCustomHbondIxn
::
AngleTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpression
.
differentiate
(
i
ter
->
first
).
optimize
().
createProgram
()));
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
dihedrals
.
begin
();
iter
!=
dihedrals
.
end
();
++
iter
)
dihedralTerms
.
push_back
(
ReferenceCustomHbondIxn
::
DihedralTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpression
.
differentiate
(
i
ter
->
first
).
optimize
().
createProgram
()));
for
(
auto
&
term
:
distances
)
distanceTerms
.
push_back
(
ReferenceCustomHbondIxn
::
DistanceTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpression
.
differentiate
(
ter
m
.
first
).
optimize
().
createProgram
()));
for
(
auto
&
term
:
angles
)
angleTerms
.
push_back
(
ReferenceCustomHbondIxn
::
AngleTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpression
.
differentiate
(
ter
m
.
first
).
optimize
().
createProgram
()));
for
(
auto
&
term
:
dihedrals
)
dihedralTerms
.
push_back
(
ReferenceCustomHbondIxn
::
DihedralTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpression
.
differentiate
(
ter
m
.
first
).
optimize
().
createProgram
()));
}
/**---------------------------------------------------------------------------------------
...
...
platforms/reference/src/SimTKReference/ReferenceCustomManyParticleIxn.cpp
View file @
08e8b206
...
...
@@ -60,8 +60,8 @@ ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(const CustomManyP
// Delete the custom functions.
for
(
map
<
string
,
Lepton
::
CustomFunction
*>::
iterator
iter
=
functions
.
begin
();
iter
!=
functions
.
end
();
iter
++
)
delete
iter
->
second
;
for
(
auto
&
function
:
functions
)
delete
function
.
second
;
// Differentiate the energy to get expressions for the force.
...
...
@@ -80,12 +80,12 @@ ReferenceCustomManyParticleIxn::ReferenceCustomManyParticleIxn(const CustomManyP
particleParamNames
[
i
].
push_back
(
paramname
.
str
());
}
}
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
distances
.
begin
();
iter
!=
distances
.
end
();
++
iter
)
distanceTerms
.
push_back
(
ReferenceCustomManyParticleIxn
::
DistanceTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpr
.
differentiate
(
i
ter
->
first
).
optimize
().
createProgram
()));
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
angles
.
begin
();
iter
!=
angles
.
end
();
++
iter
)
angleTerms
.
push_back
(
ReferenceCustomManyParticleIxn
::
AngleTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpr
.
differentiate
(
i
ter
->
first
).
optimize
().
createProgram
()));
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
dihedrals
.
begin
();
iter
!=
dihedrals
.
end
();
++
iter
)
dihedralTerms
.
push_back
(
ReferenceCustomManyParticleIxn
::
DihedralTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpr
.
differentiate
(
i
ter
->
first
).
optimize
().
createProgram
()));
for
(
auto
&
term
:
distances
)
distanceTerms
.
push_back
(
ReferenceCustomManyParticleIxn
::
DistanceTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpr
.
differentiate
(
ter
m
.
first
).
optimize
().
createProgram
()));
for
(
auto
&
term
:
angles
)
angleTerms
.
push_back
(
ReferenceCustomManyParticleIxn
::
AngleTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpr
.
differentiate
(
ter
m
.
first
).
optimize
().
createProgram
()));
for
(
auto
&
term
:
dihedrals
)
dihedralTerms
.
push_back
(
ReferenceCustomManyParticleIxn
::
DihedralTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpr
.
differentiate
(
ter
m
.
first
).
optimize
().
createProgram
()));
// Record exclusions.
...
...
@@ -192,23 +192,18 @@ void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particle
// Compute all of the variables the energy can depend on.
for
(
int
i
=
0
;
i
<
(
int
)
particleTerms
.
size
();
i
++
)
{
const
ParticleTermInfo
&
term
=
particleTerms
[
i
];
for
(
auto
&
term
:
particleTerms
)
variables
[
term
.
name
]
=
atomCoordinates
[
permutedParticles
[
term
.
atom
]][
term
.
component
];
}
for
(
int
i
=
0
;
i
<
(
int
)
distanceTerms
.
size
();
i
++
)
{
const
DistanceTermInfo
&
term
=
distanceTerms
[
i
];
for
(
auto
&
term
:
distanceTerms
)
{
computeDelta
(
permutedParticles
[
term
.
p1
],
permutedParticles
[
term
.
p2
],
term
.
delta
,
atomCoordinates
);
variables
[
term
.
name
]
=
term
.
delta
[
ReferenceForce
::
RIndex
];
}
for
(
int
i
=
0
;
i
<
(
int
)
angleTerms
.
size
();
i
++
)
{
const
AngleTermInfo
&
term
=
angleTerms
[
i
];
for
(
auto
&
term
:
angleTerms
)
{
computeDelta
(
permutedParticles
[
term
.
p1
],
permutedParticles
[
term
.
p2
],
term
.
delta1
,
atomCoordinates
);
computeDelta
(
permutedParticles
[
term
.
p3
],
permutedParticles
[
term
.
p2
],
term
.
delta2
,
atomCoordinates
);
variables
[
term
.
name
]
=
computeAngle
(
term
.
delta1
,
term
.
delta2
);
}
for
(
int
i
=
0
;
i
<
(
int
)
dihedralTerms
.
size
();
i
++
)
{
const
DihedralTermInfo
&
term
=
dihedralTerms
[
i
];
for
(
auto
&
term
:
dihedralTerms
)
{
computeDelta
(
permutedParticles
[
term
.
p2
],
permutedParticles
[
term
.
p1
],
term
.
delta1
,
atomCoordinates
);
computeDelta
(
permutedParticles
[
term
.
p2
],
permutedParticles
[
term
.
p3
],
term
.
delta2
,
atomCoordinates
);
computeDelta
(
permutedParticles
[
term
.
p4
],
permutedParticles
[
term
.
p3
],
term
.
delta3
,
atomCoordinates
);
...
...
@@ -219,15 +214,12 @@ void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particle
// Apply forces based on individual particle coordinates.
for
(
int
i
=
0
;
i
<
(
int
)
particleTerms
.
size
();
i
++
)
{
const
ParticleTermInfo
&
term
=
particleTerms
[
i
];
for
(
auto
&
term
:
particleTerms
)
forces
[
permutedParticles
[
term
.
atom
]][
term
.
component
]
-=
term
.
forceExpression
.
evaluate
(
variables
);
}
// Apply forces based on distances.
for
(
int
i
=
0
;
i
<
(
int
)
distanceTerms
.
size
();
i
++
)
{
const
DistanceTermInfo
&
term
=
distanceTerms
[
i
];
for
(
auto
&
term
:
distanceTerms
)
{
double
dEdR
=
term
.
forceExpression
.
evaluate
(
variables
)
/
(
term
.
delta
[
ReferenceForce
::
RIndex
]);
for
(
int
i
=
0
;
i
<
3
;
i
++
)
{
double
force
=
-
dEdR
*
term
.
delta
[
i
];
...
...
@@ -238,8 +230,7 @@ void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particle
// Apply forces based on angles.
for
(
int
i
=
0
;
i
<
(
int
)
angleTerms
.
size
();
i
++
)
{
const
AngleTermInfo
&
term
=
angleTerms
[
i
];
for
(
auto
&
term
:
angleTerms
)
{
double
dEdTheta
=
term
.
forceExpression
.
evaluate
(
variables
);
double
thetaCross
[
ReferenceForce
::
LastDeltaRIndex
];
SimTKOpenMMUtilities
::
crossProductVector3
(
term
.
delta1
,
term
.
delta2
,
thetaCross
);
...
...
@@ -265,8 +256,7 @@ void ReferenceCustomManyParticleIxn::calculateOneIxn(const vector<int>& particle
// Apply forces based on dihedrals.
for
(
int
i
=
0
;
i
<
(
int
)
dihedralTerms
.
size
();
i
++
)
{
const
DihedralTermInfo
&
term
=
dihedralTerms
[
i
];
for
(
auto
&
term
:
dihedralTerms
)
{
double
dEdTheta
=
term
.
forceExpression
.
evaluate
(
variables
);
double
internalF
[
4
][
3
];
double
forceFactors
[
4
];
...
...
platforms/reference/src/SimTKReference/ReferenceCustomNonbondedIxn.cpp
View file @
08e8b206
...
...
@@ -53,10 +53,10 @@ ReferenceCustomNonbondedIxn::ReferenceCustomNonbondedIxn(const Lepton::CompiledE
for
(
int
i
=
0
;
i
<
this
->
energyParamDerivExpressions
.
size
();
i
++
)
expressionSet
.
registerExpression
(
this
->
energyParamDerivExpressions
[
i
]);
rIndex
=
expressionSet
.
getVariableIndex
(
"r"
);
for
(
int
i
=
0
;
i
<
(
int
)
paramNames
.
size
();
i
++
)
{
for
(
auto
&
param
:
paramNames
)
{
for
(
int
j
=
1
;
j
<
3
;
j
++
)
{
stringstream
name
;
name
<<
param
Names
[
i
]
<<
j
;
name
<<
param
<<
j
;
particleParamIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
name
.
str
()));
}
}
...
...
@@ -159,14 +159,14 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn(int numberOfAtoms, vector<Vec
double
*
fixedParameters
,
const
map
<
string
,
double
>&
globalParameters
,
vector
<
Vec3
>&
forces
,
double
*
energyByAtom
,
double
*
totalEnergy
,
double
*
energyParamDerivs
)
{
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globalParameters
.
begin
();
iter
!=
globalParameters
.
end
();
++
iter
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
if
(
interactionGroups
.
size
()
>
0
)
{
for
(
auto
&
param
:
globalParameters
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
param
.
first
),
param
.
second
);
if
(
interactionGroups
.
size
()
>
0
)
{
// The user has specified interaction groups, so compute only the requested interactions.
for
(
int
group
=
0
;
group
<
(
int
)
interactionGroups
.
size
();
group
++
)
{
const
set
<
int
>&
set1
=
interactionGroups
[
group
]
.
first
;
const
set
<
int
>&
set2
=
interactionGroups
[
group
]
.
second
;
for
(
auto
&
group
:
interactionGroups
)
{
const
set
<
int
>&
set1
=
group
.
first
;
const
set
<
int
>&
set2
=
group
.
second
;
for
(
set
<
int
>::
const_iterator
atom1
=
set1
.
begin
();
atom1
!=
set1
.
end
();
++
atom1
)
{
for
(
set
<
int
>::
const_iterator
atom2
=
set2
.
begin
();
atom2
!=
set2
.
end
();
++
atom2
)
{
if
(
*
atom1
==
*
atom2
||
exclusions
[
*
atom1
].
find
(
*
atom2
)
!=
exclusions
[
*
atom1
].
end
())
...
...
@@ -185,8 +185,7 @@ void ReferenceCustomNonbondedIxn::calculatePairIxn(int numberOfAtoms, vector<Vec
else
if
(
cutoff
)
{
// We are using a cutoff, so get the interactions from the neighbor list.
for
(
int
i
=
0
;
i
<
(
int
)
neighborList
->
size
();
i
++
)
{
OpenMM
::
AtomPair
pair
=
(
*
neighborList
)[
i
];
for
(
auto
&
pair
:
*
neighborList
)
{
for
(
int
j
=
0
;
j
<
(
int
)
paramNames
.
size
();
j
++
)
{
expressionSet
.
setVariable
(
particleParamIndex
[
j
*
2
],
atomParameters
[
pair
.
first
][
j
]);
expressionSet
.
setVariable
(
particleParamIndex
[
j
*
2
+
1
],
atomParameters
[
pair
.
second
][
j
]);
...
...
platforms/reference/src/SimTKReference/ReferenceCustomTorsionIxn.cpp
View file @
08e8b206
...
...
@@ -47,10 +47,10 @@ ReferenceCustomTorsionIxn::ReferenceCustomTorsionIxn(const Lepton::CompiledExpre
expressionSet
.
registerExpression
(
this
->
energyParamDerivExpressions
[
i
]);
thetaIndex
=
expressionSet
.
getVariableIndex
(
"theta"
);
numParameters
=
parameterNames
.
size
();
for
(
int
i
=
0
;
i
<
(
int
)
numP
arameter
s
;
i
++
)
torsionParamIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
param
eterNames
[
i
]
));
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globalParameters
.
begin
();
iter
!=
globalParameters
.
end
();
++
iter
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
auto
&
param
:
p
arameter
Names
)
torsionParamIndex
.
push_back
(
expressionSet
.
getVariableIndex
(
param
));
for
(
auto
&
param
:
globalParameters
)
expressionSet
.
setVariable
(
expressionSet
.
getVariableIndex
(
param
.
first
),
param
.
second
);
}
/**---------------------------------------------------------------------------------------
...
...
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
View file @
08e8b206
...
...
@@ -395,8 +395,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<Vec3>& a
double
totalRealSpaceEwaldEnergy
=
0.0
f
;
for
(
int
i
=
0
;
i
<
(
int
)
neighborList
->
size
();
i
++
)
{
OpenMM
::
AtomPair
pair
=
(
*
neighborList
)[
i
];
for
(
auto
&
pair
:
*
neighborList
)
{
int
ii
=
pair
.
first
;
int
jj
=
pair
.
second
;
...
...
@@ -489,10 +488,10 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<Vec3>& a
double
totalExclusionEnergy
=
0.0
f
;
const
double
TWO_OVER_SQRT_PI
=
2
/
sqrt
(
PI_M
);
for
(
int
i
=
0
;
i
<
numberOfAtoms
;
i
++
)
for
(
set
<
int
>::
const_iterator
iter
=
exclusions
[
i
].
begin
();
iter
!=
exclusions
[
i
].
end
();
++
iter
)
{
if
(
*
iter
>
i
)
{
for
(
int
exclusion
:
exclusions
[
i
]
)
{
if
(
exclusion
>
i
)
{
int
ii
=
i
;
int
jj
=
*
iter
;
int
jj
=
exclusion
;
double
deltaR
[
2
][
ReferenceForce
::
LastDeltaRIndex
];
ReferenceForce
::
getDeltaR
(
atomCoordinates
[
jj
],
atomCoordinates
[
ii
],
deltaR
[
0
]);
...
...
@@ -581,10 +580,8 @@ void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<Vec3>& at
if
(
!
includeDirect
)
return
;
if
(
cutoff
)
{
for
(
int
i
=
0
;
i
<
(
int
)
neighborList
->
size
();
i
++
)
{
OpenMM
::
AtomPair
pair
=
(
*
neighborList
)[
i
];
for
(
auto
&
pair
:
*
neighborList
)
calculateOneIxn
(
pair
.
first
,
pair
.
second
,
atomCoordinates
,
atomParameters
,
forces
,
energyByAtom
,
totalEnergy
);
}
}
else
{
for
(
int
ii
=
0
;
ii
<
numberOfAtoms
;
ii
++
)
{
...
...
platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp
View file @
08e8b206
...
...
@@ -72,15 +72,15 @@ void ReferenceMonteCarloBarostat::applyBarostat(vector<Vec3>& atomPositions, con
// Loop over molecules.
for
(
int
i
=
0
;
i
<
(
int
)
molecules
.
size
();
i
++
)
{
for
(
auto
&
molecule
:
molecules
)
{
// Find the molecule center.
Vec3
pos
(
0
,
0
,
0
);
for
(
int
j
=
0
;
j
<
(
int
)
molecules
[
i
].
size
();
j
++
)
{
Vec3
&
atomPos
=
atomPositions
[
molecules
[
i
][
j
]
];
for
(
int
atom
:
molecule
)
{
Vec3
&
atomPos
=
atomPositions
[
atom
];
pos
+=
atomPos
;
}
pos
/=
molecule
s
[
i
]
.
size
();
pos
/=
molecule
.
size
();
// Move it into the first periodic box.
...
...
@@ -95,8 +95,8 @@ void ReferenceMonteCarloBarostat::applyBarostat(vector<Vec3>& atomPositions, con
newPos
[
1
]
*=
scaleY
;
newPos
[
2
]
*=
scaleZ
;
Vec3
offset
=
newPos
-
pos
;
for
(
int
j
=
0
;
j
<
(
int
)
molecules
[
i
].
size
();
j
++
)
{
Vec3
&
atomPos
=
atomPositions
[
molecules
[
i
][
j
]
];
for
(
int
atom
:
molecule
)
{
Vec3
&
atomPos
=
atomPositions
[
atom
];
atomPos
+=
offset
;
}
}
...
...
platforms/reference/src/SimTKReference/ReferenceNeighborList.cpp
View file @
08e8b206
...
...
@@ -184,10 +184,10 @@ public:
const
map
<
VoxelIndex
,
Voxel
>::
const_iterator
voxelEntry
=
voxelMap
.
find
(
voxelIndex
);
if
(
voxelEntry
==
voxelMap
.
end
())
continue
;
// no such voxel; skip
const
Voxel
&
voxel
=
voxelEntry
->
second
;
for
(
Voxel
::
const_iter
ato
r
item
Iter
=
voxel
.
begin
();
itemIter
!=
voxel
.
end
();
++
itemIter
)
for
(
a
u
to
&
item
:
voxel
)
{
const
AtomIndex
atomJ
=
item
Iter
->
second
;
const
Vec3
&
locationJ
=
*
item
Iter
->
first
;
const
AtomIndex
atomJ
=
item
.
second
;
const
Vec3
&
locationJ
=
*
item
.
first
;
// Ignore self hits
if
(
atomI
==
atomJ
)
continue
;
...
...
plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp
View file @
08e8b206
...
...
@@ -161,14 +161,14 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
// Double loop over different atom types.
std
::
string
sigmaCombiningRule
=
force
.
getSigmaCombiningRule
();
std
::
string
epsilonCombiningRule
=
force
.
getEpsilonCombiningRule
();
for
(
map
<
pair
<
double
,
double
>
,
int
>::
const_iterator
class1
=
classCounts
.
begin
();
class1
!=
classCounts
.
end
();
++
class1
)
{
for
(
auto
&
class1
:
classCounts
)
{
k
=
0
;
for
(
map
<
pair
<
double
,
double
>
,
int
>::
const_iterator
class2
=
classCounts
.
begin
();
class2
!=
classCounts
.
end
();
++
class2
)
{
for
(
auto
&
class2
:
classCounts
)
{
// AMOEBA combining rules, copied over from the CUDA code.
double
iSigma
=
class1
->
first
.
first
;
double
jSigma
=
class2
->
first
.
first
;
double
iEpsilon
=
class1
->
first
.
second
;
double
jEpsilon
=
class2
->
first
.
second
;
double
iSigma
=
class1
.
first
.
first
;
double
jSigma
=
class2
.
first
.
first
;
double
iEpsilon
=
class1
.
first
.
second
;
double
jEpsilon
=
class2
.
first
.
second
;
// ARITHMETIC = 1
// GEOMETRIC = 2
// CUBIC-MEAN = 3
...
...
@@ -207,7 +207,7 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
epsilon
=
0.0
;
}
}
int
count
=
class1
->
second
*
class2
->
second
;
int
count
=
class1
.
second
*
class2
.
second
;
// Below is an exact copy of stuff from the previous block.
double
rv
=
sigma
;
double
termik
=
2.0
*
M_PI
*
count
;
// termik is equivalent to 2 * pi * count.
...
...
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @
08e8b206
...
...
@@ -973,8 +973,8 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
molecularQuadrupolesVec
.
push_back
((
float
)
quadrupole
[
5
]);
}
hasQuadrupoles
=
false
;
for
(
int
i
=
0
;
i
<
(
int
)
molecularQuadrupolesVec
.
size
();
i
++
)
if
(
molecularQuadrupolesVec
[
i
]
!=
0.0
)
for
(
auto
q
:
molecularQuadrupolesVec
)
if
(
q
!=
0.0
)
hasQuadrupoles
=
true
;
int
paddedNumAtoms
=
cu
.
getPaddedNumAtoms
();
for
(
int
i
=
numMultipoles
;
i
<
paddedNumAtoms
;
i
++
)
{
...
...
@@ -1049,15 +1049,15 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
Covalent13
,
atoms
);
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
for
(
set
<
int
>::
const_iterator
iter
=
allAtoms
.
begin
();
iter
!=
allAtoms
.
end
();
++
iter
)
covalentFlagValues
.
push_back
(
make_int3
(
i
,
*
iter
,
0
));
for
(
int
atom
:
allAtoms
)
covalentFlagValues
.
push_back
(
make_int3
(
i
,
atom
,
0
));
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
Covalent14
,
atoms
);
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
for
(
int
j
=
0
;
j
<
(
int
)
atoms
.
size
();
j
++
)
covalentFlagValues
.
push_back
(
make_int3
(
i
,
atom
s
[
j
]
,
1
));
for
(
int
atom
:
atoms
)
covalentFlagValues
.
push_back
(
make_int3
(
i
,
atom
,
1
));
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
Covalent15
,
atoms
);
for
(
int
j
=
0
;
j
<
(
int
)
atoms
.
size
();
j
++
)
covalentFlagValues
.
push_back
(
make_int3
(
i
,
atom
s
[
j
]
,
2
));
for
(
int
atom
:
atoms
)
covalentFlagValues
.
push_back
(
make_int3
(
i
,
atom
,
2
));
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
PolarizationCovalent11
,
atoms
);
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
...
...
@@ -1068,15 +1068,14 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
vector
<
int
>
atoms12
;
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
PolarizationCovalent12
,
atoms12
);
for
(
int
j
=
0
;
j
<
(
int
)
atoms
.
size
();
j
++
)
if
(
find
(
atoms12
.
begin
(),
atoms12
.
end
(),
atom
s
[
j
]
)
==
atoms12
.
end
())
polarizationFlagValues
.
push_back
(
make_int2
(
i
,
atom
s
[
j
]
));
for
(
int
atom
:
atoms
)
if
(
find
(
atoms12
.
begin
(),
atoms12
.
end
(),
atom
)
==
atoms12
.
end
())
polarizationFlagValues
.
push_back
(
make_int2
(
i
,
atom
));
}
set
<
pair
<
int
,
int
>
>
tilesWithExclusions
;
for
(
int
atom1
=
0
;
atom1
<
(
int
)
exclusions
.
size
();
++
atom1
)
{
int
x
=
atom1
/
CudaContext
::
TileSize
;
for
(
int
j
=
0
;
j
<
(
int
)
exclusions
[
atom1
].
size
();
++
j
)
{
int
atom2
=
exclusions
[
atom1
][
j
];
for
(
int
atom2
:
exclusions
[
atom1
])
{
int
y
=
atom2
/
CudaContext
::
TileSize
;
tilesWithExclusions
.
insert
(
make_pair
(
max
(
x
,
y
),
min
(
x
,
y
)));
}
...
...
@@ -1412,10 +1411,10 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
}
covalentFlags
=
CudaArray
::
create
<
uint2
>
(
cu
,
nb
.
getExclusions
().
getSize
(),
"covalentFlags"
);
vector
<
uint2
>
covalentFlagsVec
(
nb
.
getExclusions
().
getSize
(),
make_uint2
(
0
,
0
));
for
(
int
i
=
0
;
i
<
(
int
)
covalentFlagValues
.
size
();
i
++
)
{
int
atom1
=
co
val
entFlagValues
[
i
]
.
x
;
int
atom2
=
co
val
entFlagValues
[
i
]
.
y
;
int
value
=
co
val
entFlagValues
[
i
]
.
z
;
for
(
int
3
values
:
covalentFlagValues
)
{
int
atom1
=
val
ues
.
x
;
int
atom2
=
val
ues
.
y
;
int
value
=
val
ues
.
z
;
int
x
=
atom1
/
CudaContext
::
TileSize
;
int
offset1
=
atom1
-
x
*
CudaContext
::
TileSize
;
int
y
=
atom2
/
CudaContext
::
TileSize
;
...
...
@@ -1446,9 +1445,9 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
polarizationGroupFlags
=
CudaArray
::
create
<
unsigned
int
>
(
cu
,
nb
.
getExclusions
().
getSize
(),
"polarizationGroupFlags"
);
vector
<
unsigned
int
>
polarizationGroupFlagsVec
(
nb
.
getExclusions
().
getSize
(),
0
);
for
(
int
i
=
0
;
i
<
(
int
)
polarizationFlagValues
.
size
();
i
++
)
{
int
atom1
=
polarizationFlagV
alues
[
i
]
.
x
;
int
atom2
=
polarizationFlagV
alues
[
i
]
.
y
;
for
(
int
2
values
:
polarizationFlagValues
)
{
int
atom1
=
v
alues
.
x
;
int
atom2
=
v
alues
.
y
;
int
x
=
atom1
/
CudaContext
::
TileSize
;
int
offset1
=
atom1
-
x
*
CudaContext
::
TileSize
;
int
y
=
atom2
/
CudaContext
::
TileSize
;
...
...
@@ -1473,10 +1472,12 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
double
CudaCalcAmoebaMultipoleForceKernel
::
execute
(
ContextImpl
&
context
,
bool
includeForces
,
bool
includeEnergy
)
{
if
(
!
hasInitializedScaleFactors
)
{
initializeScaleFactors
();
for
(
int
i
=
0
;
i
<
(
int
)
context
.
getForceImpls
()
.
size
()
&&
gkKernel
==
NULL
;
i
++
)
{
AmoebaGeneralizedKirkwoodForceImpl
*
gkImpl
=
dynamic_cast
<
AmoebaGeneralizedKirkwoodForceImpl
*>
(
context
.
getForceImpls
()[
i
]
);
if
(
gkImpl
!=
NULL
)
for
(
auto
impl
:
context
.
getForceImpls
())
{
AmoebaGeneralizedKirkwoodForceImpl
*
gkImpl
=
dynamic_cast
<
AmoebaGeneralizedKirkwoodForceImpl
*>
(
impl
);
if
(
gkImpl
!=
NULL
)
{
gkKernel
=
dynamic_cast
<
CudaCalcAmoebaGeneralizedKirkwoodForceKernel
*>
(
&
gkImpl
->
getKernel
().
getImpl
());
break
;
}
}
}
CudaNonbondedUtilities
&
nb
=
cu
.
getNonbondedUtilities
();
...
...
@@ -2232,8 +2233,8 @@ void CudaCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImpl& co
molecularQuadrupolesVec
.
push_back
((
float
)
quadrupole
[
5
]);
}
if
(
!
hasQuadrupoles
)
{
for
(
int
i
=
0
;
i
<
(
int
)
molecularQuadrupolesVec
.
size
();
i
++
)
if
(
molecularQuadrupolesVec
[
i
]
!=
0.0
)
for
(
auto
q
:
molecularQuadrupolesVec
)
if
(
q
!=
0.0
)
throw
OpenMMException
(
"updateParametersInContext: Cannot set a non-zero quadrupole moment, because quadrupoles were excluded from the kernel"
);
}
for
(
int
i
=
force
.
getNumMultipoles
();
i
<
cu
.
getPaddedNumAtoms
();
i
++
)
{
...
...
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaExtrapolatedPolarization.cpp
View file @
08e8b206
...
...
@@ -290,8 +290,8 @@ static void check_finite_differences(vector<Vec3> analytic_forces, Context &cont
// 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
<
(
int
)
analytic_forces
.
size
();
++
i
)
norm
+=
analytic_forces
[
i
].
dot
(
analytic_forces
[
i
]
);
for
(
auto
&
f
:
analytic_forces
)
norm
+=
f
.
dot
(
f
);
norm
=
std
::
sqrt
(
norm
);
const
double
stepSize
=
1e-3
;
double
step
=
0.5
*
stepSize
/
norm
;
...
...
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp
View file @
08e8b206
...
...
@@ -61,8 +61,8 @@ static void checkFiniteDifferences(vector<Vec3> forces, Context &context, vector
// 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 < (int) forces.size(); ++i
)
norm += f
orces[i].dot(forces[i]
);
for (
auto& f : forces
)
norm += f
.dot(f
);
norm = std::sqrt(norm);
const double stepSize = 1e-3;
double step = 0.5*stepSize/norm;
...
...
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaMultipoleForce.cpp
View file @
08e8b206
...
...
@@ -2919,8 +2919,8 @@ static void testNoQuadrupoles(bool usePme) {
int
axisType
,
atomX
,
atomY
,
atomZ
;
vector
<
double
>
dipole
,
quadrupole
;
amoebaMultipoleForce
->
getMultipoleParameters
(
i
,
charge
,
dipole
,
quadrupole
,
axisType
,
atomZ
,
atomX
,
atomY
,
thole
,
damping
,
polarity
);
for
(
int
j
=
0
;
j
<
(
int
)
quadrupole
.
size
();
j
++
)
q
uadrupole
[
j
]
=
0
;
for
(
auto
&
q
:
quadrupole
)
q
=
0
;
amoebaMultipoleForce
->
setMultipoleParameters
(
i
,
charge
,
dipole
,
quadrupole
,
axisType
,
atomZ
,
atomX
,
atomY
,
thole
,
damping
,
polarity
);
}
amoebaMultipoleForce
->
updateParametersInContext
(
context
);
...
...
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp
View file @
08e8b206
This diff is collapsed.
Click to expand it.
Prev
1
2
3
4
5
Next
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment