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
28b79b2f
Commit
28b79b2f
authored
Mar 03, 2017
by
peastman
Browse files
Use C++11 style loops
parent
ab8f1021
Changes
21
Show whitespace changes
Inline
Side-by-side
Showing
20 changed files
with
211 additions
and
261 deletions
+211
-261
platforms/cpu/src/CpuCustomGBForce.cpp
platforms/cpu/src/CpuCustomGBForce.cpp
+50
-50
platforms/cpu/src/CpuCustomManyParticleForce.cpp
platforms/cpu/src/CpuCustomManyParticleForce.cpp
+26
-35
platforms/cpu/src/CpuCustomNonbondedForce.cpp
platforms/cpu/src/CpuCustomNonbondedForce.cpp
+12
-12
platforms/cpu/src/CpuRandom.cpp
platforms/cpu/src/CpuRandom.cpp
+2
-2
platforms/cpu/src/CpuSETTLE.cpp
platforms/cpu/src/CpuSETTLE.cpp
+2
-2
platforms/reference/src/SimTKReference/ReferenceAndersenThermostat.cpp
...erence/src/SimTKReference/ReferenceAndersenThermostat.cpp
+2
-3
platforms/reference/src/SimTKReference/ReferenceCCMAAlgorithm.cpp
...s/reference/src/SimTKReference/ReferenceCCMAAlgorithm.cpp
+2
-5
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
+2
-5
platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp
...erence/src/SimTKReference/ReferenceMonteCarloBarostat.cpp
+6
-6
No files found.
platforms/cpu/src/CpuCustomGBForce.cpp
View file @
28b79b2f
...
...
@@ -74,48 +74,48 @@ CpuCustomGBForce::ThreadData::ThreadData(int numAtoms, int numThreads, int threa
variableLocations
[
name
.
str
()]
=
&
particleValue
[
2
*
i
+
j
];
}
}
for
(
int
i
=
0
;
i
<
(
int
)
valueExpressions
.
size
();
i
++
)
{
this
->
valueE
xpression
s
[
i
]
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
this
->
valueE
xpression
s
[
i
]
);
}
for
(
int
i
=
0
;
i
<
(
int
)
valueDerivExpressions
.
size
();
i
++
)
for
(
int
j
=
0
;
j
<
(
int
)
valueDerivExpressions
[
i
].
size
();
j
++
)
{
this
->
valueDerivE
xpression
s
[
i
][
j
]
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
this
->
valueDerivE
xpression
s
[
i
][
j
]
);
}
for
(
int
i
=
0
;
i
<
(
int
)
valueGradientExpressions
.
size
();
i
++
)
for
(
int
j
=
0
;
j
<
(
int
)
valueGradientExpressions
[
i
].
size
();
j
++
)
{
this
->
valueGradientE
xpression
s
[
i
][
j
]
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
this
->
valueGradientE
xpression
s
[
i
][
j
]
);
}
for
(
int
i
=
0
;
i
<
(
int
)
valueParamDerivExpressions
.
size
();
i
++
)
for
(
int
j
=
0
;
j
<
(
int
)
valueParamDerivExpressions
[
i
].
size
();
j
++
)
{
this
->
valueParamDerivE
xpression
s
[
i
][
j
]
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
this
->
valueParamDerivE
xpression
s
[
i
][
j
]
);
}
for
(
int
i
=
0
;
i
<
(
int
)
energyExpressions
.
size
();
i
++
)
{
this
->
energyE
xpression
s
[
i
]
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
this
->
energyE
xpression
s
[
i
]
);
}
for
(
int
i
=
0
;
i
<
(
int
)
energyDerivExpressions
.
size
();
i
++
)
for
(
int
j
=
0
;
j
<
(
int
)
energyDerivExpressions
[
i
].
size
();
j
++
)
{
this
->
energyDerivE
xpression
s
[
i
][
j
]
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
this
->
energyDerivE
xpression
s
[
i
][
j
]
);
}
for
(
int
i
=
0
;
i
<
(
int
)
energyGradientExpressions
.
size
();
i
++
)
for
(
int
j
=
0
;
j
<
(
int
)
energyGradientExpressions
[
i
].
size
();
j
++
)
{
this
->
energyGradientE
xpression
s
[
i
][
j
]
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
this
->
energyGradientE
xpression
s
[
i
][
j
]
);
}
for
(
int
i
=
0
;
i
<
(
int
)
energyParamDerivExpressions
.
size
();
i
++
)
for
(
int
j
=
0
;
j
<
(
int
)
energyParamDerivExpressions
[
i
].
size
();
j
++
)
{
this
->
energyParamDerivE
xpression
s
[
i
][
j
]
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
this
->
energyParamDerivE
xpression
s
[
i
][
j
]
);
for
(
auto
&
expression
:
this
->
valueExpressions
)
{
e
xpression
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
e
xpression
);
}
for
(
auto
&
expressions
:
this
->
valueDerivExpressions
)
for
(
auto
&
expression
:
expressions
)
{
e
xpression
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
e
xpression
);
}
for
(
auto
&
expressions
:
this
->
valueGradientExpressions
)
for
(
auto
&
expression
:
expressions
)
{
e
xpression
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
e
xpression
);
}
for
(
auto
&
expressions
:
this
->
valueParamDerivExpressions
)
for
(
auto
&
expression
:
expressions
)
{
e
xpression
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
e
xpression
);
}
for
(
auto
&
expression
:
this
->
energyExpressions
)
{
e
xpression
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
e
xpression
);
}
for
(
auto
&
expressions
:
this
->
energyDerivExpressions
)
for
(
auto
&
expression
:
expressions
)
{
e
xpression
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
e
xpression
);
}
for
(
auto
&
expressions
:
this
->
energyGradientExpressions
)
for
(
auto
&
expression
:
expressions
)
{
e
xpression
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
e
xpression
);
}
for
(
auto
&
expressions
:
this
->
energyParamDerivExpressions
)
for
(
auto
&
expression
:
expressions
)
{
e
xpression
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
e
xpression
);
}
value0
.
resize
(
numAtoms
);
dEdV
.
resize
(
valueNames
.
size
());
for
(
int
i
=
0
;
i
<
(
int
)
dEdV
.
size
();
i
++
)
dEdV
[
i
]
.
resize
(
numAtoms
);
for
(
auto
&
v
:
dEdV
)
v
.
resize
(
numAtoms
);
dVdX
.
resize
(
valueDerivExpressions
.
size
());
dVdY
.
resize
(
valueDerivExpressions
.
size
());
dVdZ
.
resize
(
valueDerivExpressions
.
size
());
...
...
@@ -155,8 +155,8 @@ CpuCustomGBForce::CpuCustomGBForce(int numAtoms, const std::vector<std::set<int>
}
CpuCustomGBForce
::~
CpuCustomGBForce
()
{
for
(
int
i
=
0
;
i
<
(
int
)
threadData
.
size
();
i
++
)
delete
threadData
[
i
]
;
for
(
auto
data
:
threadData
)
delete
data
;
}
void
CpuCustomGBForce
::
setUseCutoff
(
float
distance
,
const
CpuNeighborList
&
neighbors
)
{
...
...
@@ -256,16 +256,16 @@ void CpuCustomGBForce::threadComputeForce(ThreadPool& threads, int threadIndex)
ThreadData
&
data
=
*
threadData
[
threadIndex
];
fvec4
boxSize
(
periodicBoxSize
[
0
],
periodicBoxSize
[
1
],
periodicBoxSize
[
2
],
0
);
fvec4
invBoxSize
((
1
/
periodicBoxSize
[
0
]),
(
1
/
periodicBoxSize
[
1
]),
(
1
/
periodicBoxSize
[
2
]),
0
);
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globalParameters
->
begin
();
iter
!=
globalParameters
->
end
();
++
iter
)
data
.
expressionSet
.
setVariable
(
data
.
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
auto
&
param
:
*
globalParameters
)
data
.
expressionSet
.
setVariable
(
data
.
expressionSet
.
getVariableIndex
(
param
.
first
),
param
.
second
);
// Calculate the first computed value.
for
(
int
i
=
0
;
i
<
(
int
)
data
.
value0
.
size
();
i
++
)
data
.
value0
[
i
]
=
0.0
f
;
for
(
int
i
=
0
;
i
<
(
int
)
data
.
dValue0dParam
.
size
();
i
++
)
for
(
int
j
=
0
;
j
<
(
int
)
data
.
dValue0dParam
[
i
].
size
();
j
++
)
data
.
dValue0dParam
[
i
][
j
]
=
0.0
;
for
(
auto
&
v
:
data
.
value0
)
v
=
0.0
f
;
for
(
auto
&
vals
:
data
.
dValue0dParam
)
for
(
auto
&
v
:
vals
)
v
=
0.0
f
;
if
(
valueTypes
[
0
]
==
CustomGBForce
::
ParticlePair
)
calculateParticlePairValue
(
0
,
data
,
numberOfAtoms
,
posq
,
atomParameters
,
true
,
boxSize
,
invBoxSize
);
else
...
...
@@ -291,8 +291,8 @@ void CpuCustomGBForce::threadComputeForce(ThreadPool& threads, int threadIndex)
int
numValues
=
valueTypes
.
size
();
for
(
int
atom
=
data
.
firstAtom
;
atom
<
data
.
lastAtom
;
atom
++
)
{
float
sum
=
0.0
f
;
for
(
int
j
=
0
;
j
<
(
int
)
threadData
.
size
();
j
++
)
sum
+=
threadData
[
j
]
->
value0
[
atom
];
for
(
auto
&
data
:
threadData
)
sum
+=
data
->
value0
[
atom
];
values
[
0
][
atom
]
=
sum
;
data
.
x
=
posq
[
4
*
atom
];
data
.
y
=
posq
[
4
*
atom
+
1
];
...
...
@@ -320,11 +320,11 @@ void CpuCustomGBForce::threadComputeForce(ThreadPool& threads, int threadIndex)
// Now calculate the energy and its derivatives.
for
(
int
i
=
0
;
i
<
(
int
)
data
.
dEdV
.
size
();
i
++
)
for
(
int
j
=
0
;
j
<
(
int
)
data
.
dEdV
[
i
].
size
();
j
++
)
data
.
dEdV
[
i
][
j
]
=
0.0
f
;
for
(
int
i
=
0
;
i
<
(
int
)
data
.
energyParamDerivs
.
size
();
i
++
)
data
.
energyParamDerivs
[
i
]
=
0.0
f
;
for
(
auto
&
vals
:
data
.
dEdV
)
for
(
auto
&
v
:
vals
)
v
=
0.0
f
;
for
(
auto
&
v
:
data
.
energyParamDerivs
)
v
=
0.0
f
;
for
(
int
termIndex
=
0
;
termIndex
<
(
int
)
data
.
energyExpressions
.
size
();
termIndex
++
)
{
if
(
energyTypes
[
termIndex
]
==
CustomGBForce
::
SingleParticle
)
calculateSingleParticleEnergyTerm
(
termIndex
,
data
,
numberOfAtoms
,
posq
,
atomParameters
,
forces
,
energy
);
...
...
@@ -340,8 +340,8 @@ void CpuCustomGBForce::threadComputeForce(ThreadPool& threads, int threadIndex)
for
(
int
atom
=
data
.
firstAtom
;
atom
<
data
.
lastAtom
;
atom
++
)
{
for
(
int
i
=
0
;
i
<
(
int
)
dEdV
.
size
();
i
++
)
{
float
sum
=
0.0
f
;
for
(
int
j
=
0
;
j
<
(
int
)
threadData
.
size
();
j
++
)
sum
+=
threadData
[
j
]
->
dEdV
[
i
][
atom
];
for
(
auto
&
data
:
threadData
)
sum
+=
data
->
dEdV
[
i
][
atom
];
dEdV
[
i
][
atom
]
=
sum
;
}
}
...
...
platforms/cpu/src/CpuCustomManyParticleForce.cpp
View file @
28b79b2f
...
...
@@ -63,8 +63,8 @@ CpuCustomManyParticleForce::CpuCustomManyParticleForce(const CustomManyParticleF
// 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
;
// Record exclusions.
...
...
@@ -84,8 +84,8 @@ CpuCustomManyParticleForce::CpuCustomManyParticleForce(const CustomManyParticleF
CpuCustomManyParticleForce
::~
CpuCustomManyParticleForce
()
{
if
(
neighborList
!=
NULL
)
delete
neighborList
;
for
(
int
i
=
0
;
i
<
(
int
)
threadData
.
size
();
i
++
)
delete
threadData
[
i
]
;
for
(
auto
data
:
threadData
)
delete
data
;
}
void
CpuCustomManyParticleForce
::
calculateIxn
(
AlignedArray
<
float
>&
posq
,
double
**
particleParameters
,
...
...
@@ -150,8 +150,8 @@ void CpuCustomManyParticleForce::threadComputeForce(ThreadPool& threads, int thr
float
*
forces
=
&
(
*
threadForce
)[
threadIndex
][
0
];
ThreadData
&
data
=
*
threadData
[
threadIndex
];
data
.
energy
=
0
;
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globalParameters
->
begin
();
iter
!=
globalParameters
->
end
();
++
iter
)
data
.
expressionSet
.
setVariable
(
data
.
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
auto
&
param
:
*
globalParameters
)
data
.
expressionSet
.
setVariable
(
data
.
expressionSet
.
getVariableIndex
(
param
.
first
),
param
.
second
);
if
(
useCutoff
)
{
// Loop over interactions from the neighbor list.
...
...
@@ -287,18 +287,12 @@ void CpuCustomManyParticleForce::calculateOneIxn(vector<int>& particleSet, doubl
// Compute all of the variables the energy can depend on.
for
(
int
i
=
0
;
i
<
(
int
)
data
.
particleTerms
.
size
();
i
++
)
{
const
ParticleTermInfo
&
term
=
data
.
particleTerms
[
i
];
for
(
auto
&
term
:
data
.
particleTerms
)
expressionSet
.
setVariable
(
term
.
variableIndex
,
posq
[
4
*
permutedParticles
[
term
.
atom
]
+
term
.
component
]);
}
for
(
int
i
=
0
;
i
<
(
int
)
data
.
distanceTerms
.
size
();
i
++
)
{
const
DistanceTermInfo
&
term
=
data
.
distanceTerms
[
i
];
for
(
auto
&
term
:
data
.
distanceTerms
)
expressionSet
.
setVariable
(
term
.
variableIndex
,
normDelta
[
term
.
delta
]);
}
for
(
int
i
=
0
;
i
<
(
int
)
data
.
angleTerms
.
size
();
i
++
)
{
const
AngleTermInfo
&
term
=
data
.
angleTerms
[
i
];
for
(
auto
&
term
:
data
.
angleTerms
)
expressionSet
.
setVariable
(
term
.
variableIndex
,
computeAngle
(
delta
[
term
.
delta1
],
delta
[
term
.
delta2
],
norm2Delta
[
term
.
delta1
],
norm2Delta
[
term
.
delta2
],
term
.
delta1Sign
*
term
.
delta2Sign
));
}
for
(
int
i
=
0
;
i
<
(
int
)
data
.
dihedralTerms
.
size
();
i
++
)
{
const
DihedralTermInfo
&
term
=
data
.
dihedralTerms
[
i
];
expressionSet
.
setVariable
(
term
.
variableIndex
,
getDihedralAngleBetweenThreeVectors
(
delta
[
term
.
delta1
],
delta
[
term
.
delta2
],
delta
[
term
.
delta3
],
cross1
[
i
],
cross2
[
i
],
delta
[
term
.
delta1
]));
...
...
@@ -310,8 +304,7 @@ void CpuCustomManyParticleForce::calculateOneIxn(vector<int>& particleSet, doubl
AlignedArray
<
fvec4
>&
f
=
data
.
f
;
for
(
int
i
=
0
;
i
<
numParticlesPerSet
;
i
++
)
f
[
i
]
=
fvec4
(
0.0
f
);
for
(
int
i
=
0
;
i
<
(
int
)
data
.
particleTerms
.
size
();
i
++
)
{
const
ParticleTermInfo
&
term
=
data
.
particleTerms
[
i
];
for
(
auto
&
term
:
data
.
particleTerms
)
{
float
temp
[
4
];
f
[
term
.
atom
].
store
(
temp
);
temp
[
term
.
component
]
-=
term
.
forceExpression
.
evaluate
();
...
...
@@ -320,8 +313,7 @@ void CpuCustomManyParticleForce::calculateOneIxn(vector<int>& particleSet, doubl
// Apply forces based on distances.
for
(
int
i
=
0
;
i
<
(
int
)
data
.
distanceTerms
.
size
();
i
++
)
{
const
DistanceTermInfo
&
term
=
data
.
distanceTerms
[
i
];
for
(
auto
&
term
:
data
.
distanceTerms
)
{
float
dEdR
=
(
float
)
(
term
.
forceExpression
.
evaluate
()
*
term
.
deltaSign
/
(
normDelta
[
term
.
delta
]));
fvec4
force
=
-
dEdR
*
delta
[
term
.
delta
];
f
[
term
.
p1
]
-=
force
;
...
...
@@ -330,8 +322,7 @@ void CpuCustomManyParticleForce::calculateOneIxn(vector<int>& particleSet, doubl
// Apply forces based on angles.
for
(
int
i
=
0
;
i
<
(
int
)
data
.
angleTerms
.
size
();
i
++
)
{
const
AngleTermInfo
&
term
=
data
.
angleTerms
[
i
];
for
(
auto
&
term
:
data
.
angleTerms
)
{
float
dEdTheta
=
(
float
)
term
.
forceExpression
.
evaluate
();
fvec4
thetaCross
=
cross
(
delta
[
term
.
delta1
],
delta
[
term
.
delta2
]);
float
lengthThetaCross
=
sqrtf
(
dot3
(
thetaCross
,
thetaCross
));
...
...
@@ -482,20 +473,20 @@ CpuCustomManyParticleForce::ThreadData::ThreadData(const CustomManyParticleForce
particleParamIndices
[
i
].
push_back
(
expressionSet
.
getVariableIndex
(
paramname
.
str
()));
}
}
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
dihedrals
.
begin
();
iter
!=
dihedrals
.
end
();
++
iter
)
dihedralTerms
.
push_back
(
CpuCustomManyParticleForce
::
DihedralTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpr
.
differentiate
(
i
ter
->
first
).
optimize
().
createCompiledExpression
(),
*
this
));
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
distances
.
begin
();
iter
!=
distances
.
end
();
++
iter
)
distanceTerms
.
push_back
(
CpuCustomManyParticleForce
::
DistanceTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpr
.
differentiate
(
i
ter
->
first
).
optimize
().
createCompiledExpression
(),
*
this
));
for
(
map
<
string
,
vector
<
int
>
>::
const_iterator
iter
=
angles
.
begin
();
iter
!=
angles
.
end
();
++
iter
)
angleTerms
.
push_back
(
CpuCustomManyParticleForce
::
AngleTermInfo
(
i
ter
->
first
,
i
ter
->
second
,
energyExpr
.
differentiate
(
i
ter
->
first
).
optimize
().
createCompiledExpression
(),
*
this
));
for
(
int
i
=
0
;
i
<
particleTerms
.
size
();
i
++
)
expressionSet
.
registerExpression
(
particleTerms
[
i
]
.
forceExpression
);
for
(
int
i
=
0
;
i
<
distanceTerms
.
size
();
i
++
)
expressionSet
.
registerExpression
(
distanceTerms
[
i
]
.
forceExpression
);
for
(
int
i
=
0
;
i
<
angleTerms
.
size
();
i
++
)
expressionSet
.
registerExpression
(
angleTerms
[
i
]
.
forceExpression
);
for
(
int
i
=
0
;
i
<
dihedralTerms
.
size
();
i
++
)
expressionSet
.
registerExpression
(
dihedralTerms
[
i
]
.
forceExpression
);
for
(
auto
&
term
:
dihedrals
)
dihedralTerms
.
push_back
(
CpuCustomManyParticleForce
::
DihedralTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpr
.
differentiate
(
ter
m
.
first
).
optimize
().
createCompiledExpression
(),
*
this
));
for
(
auto
&
term
:
distances
)
distanceTerms
.
push_back
(
CpuCustomManyParticleForce
::
DistanceTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpr
.
differentiate
(
ter
m
.
first
).
optimize
().
createCompiledExpression
(),
*
this
));
for
(
auto
&
term
:
angles
)
angleTerms
.
push_back
(
CpuCustomManyParticleForce
::
AngleTermInfo
(
ter
m
.
first
,
ter
m
.
second
,
energyExpr
.
differentiate
(
ter
m
.
first
).
optimize
().
createCompiledExpression
(),
*
this
));
for
(
auto
&
term
:
particleTerms
)
expressionSet
.
registerExpression
(
term
.
forceExpression
);
for
(
auto
&
term
:
distanceTerms
)
expressionSet
.
registerExpression
(
term
.
forceExpression
);
for
(
auto
&
term
:
angleTerms
)
expressionSet
.
registerExpression
(
term
.
forceExpression
);
for
(
auto
&
term
:
dihedralTerms
)
expressionSet
.
registerExpression
(
term
.
forceExpression
);
int
numDeltas
=
deltaPairs
.
size
();
delta
.
resize
(
numDeltas
);
normDelta
.
resize
(
numDeltas
);
...
...
platforms/cpu/src/CpuCustomNonbondedForce.cpp
View file @
28b79b2f
...
...
@@ -51,9 +51,9 @@ CpuCustomNonbondedForce::ThreadData::ThreadData(const Lepton::CompiledExpression
this
->
forceExpression
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
this
->
energyExpression
);
expressionSet
.
registerExpression
(
this
->
forceExpression
);
for
(
int
i
=
0
;
i
<
this
->
energyParamDerivExpressions
.
size
();
i
++
)
{
this
->
energyParamDerivE
xpression
s
[
i
]
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
this
->
energyParamDerivE
xpression
s
[
i
]
);
for
(
auto
&
expression
:
this
->
energyParamDerivExpressions
)
{
e
xpression
.
setVariableLocations
(
variableLocations
);
expressionSet
.
registerExpression
(
e
xpression
);
}
}
...
...
@@ -66,8 +66,8 @@ CpuCustomNonbondedForce::CpuCustomNonbondedForce(const Lepton::CompiledExpressio
}
CpuCustomNonbondedForce
::~
CpuCustomNonbondedForce
()
{
for
(
int
i
=
0
;
i
<
(
int
)
threadData
.
size
();
i
++
)
delete
threadData
[
i
]
;
for
(
auto
data
:
threadData
)
delete
data
;
}
void
CpuCustomNonbondedForce
::
setUseCutoff
(
double
distance
,
const
CpuNeighborList
&
neighbors
)
{
...
...
@@ -78,9 +78,9 @@ void CpuCustomNonbondedForce::setUseCutoff(double distance, const CpuNeighborLis
void
CpuCustomNonbondedForce
::
setInteractionGroups
(
const
vector
<
pair
<
set
<
int
>
,
set
<
int
>
>
>&
groups
)
{
useInteractionGroups
=
true
;
for
(
int
group
=
0
;
group
<
(
int
)
groups
.
size
();
group
++
)
{
const
set
<
int
>&
set1
=
group
s
[
group
]
.
first
;
const
set
<
int
>&
set2
=
group
s
[
group
]
.
second
;
for
(
auto
&
group
:
groups
)
{
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
())
...
...
@@ -167,10 +167,10 @@ void CpuCustomNonbondedForce::threadComputeForce(ThreadPool& threads, int thread
double
&
energy
=
threadEnergy
[
threadIndex
];
float
*
forces
=
&
(
*
threadForce
)[
threadIndex
][
0
];
ThreadData
&
data
=
*
threadData
[
threadIndex
];
for
(
map
<
string
,
double
>::
const_iterator
iter
=
globalParameters
->
begin
();
iter
!=
globalParameters
->
end
();
++
iter
)
data
.
expressionSet
.
setVariable
(
data
.
expressionSet
.
getVariableIndex
(
iter
->
first
),
iter
->
second
);
for
(
int
i
=
0
;
i
<
data
.
energyParamDerivs
.
size
();
i
++
)
d
ata
.
energyParamDerivs
[
i
]
=
0.0
;
for
(
auto
&
param
:
*
globalParameters
)
data
.
expressionSet
.
setVariable
(
data
.
expressionSet
.
getVariableIndex
(
param
.
first
),
param
.
second
);
for
(
auto
&
deriv
:
data
.
energyParamDerivs
)
d
eriv
=
0.0
;
fvec4
boxSize
(
periodicBoxVectors
[
0
][
0
],
periodicBoxVectors
[
1
][
1
],
periodicBoxVectors
[
2
][
2
],
0
);
fvec4
invBoxSize
(
recipBoxSize
[
0
],
recipBoxSize
[
1
],
recipBoxSize
[
2
],
0
);
if
(
useInteractionGroups
)
{
...
...
platforms/cpu/src/CpuRandom.cpp
View file @
28b79b2f
...
...
@@ -34,8 +34,8 @@ CpuRandom::CpuRandom() : hasInitialized(false) {
}
CpuRandom
::~
CpuRandom
()
{
for
(
int
i
=
0
;
i
<
(
int
)
threadRandom
.
size
();
i
++
)
delete
threadR
andom
[
i
]
;
for
(
auto
random
:
threadRandom
)
delete
r
andom
;
}
void
CpuRandom
::
initialize
(
int
seed
,
int
numThreads
)
{
...
...
platforms/cpu/src/CpuSETTLE.cpp
View file @
28b79b2f
...
...
@@ -56,8 +56,8 @@ CpuSETTLE::CpuSETTLE(const System& system, const ReferenceSETTLEAlgorithm& settl
}
CpuSETTLE
::~
CpuSETTLE
()
{
for
(
int
i
=
0
;
i
<
(
int
)
threadSettle
.
size
();
i
++
)
delete
threadS
ettle
[
i
]
;
for
(
auto
settle
:
threadSettle
)
delete
s
ettle
;
}
void
CpuSETTLE
::
apply
(
vector
<
OpenMM
::
Vec3
>&
atomCoordinates
,
vector
<
OpenMM
::
Vec3
>&
atomCoordinatesP
,
vector
<
double
>&
inverseMasses
,
double
tolerance
)
{
...
...
platforms/reference/src/SimTKReference/ReferenceAndersenThermostat.cpp
View file @
28b79b2f
...
...
@@ -66,13 +66,12 @@ using namespace OpenMM;
double
temperature
,
double
collisionFrequency
,
double
stepSize
)
const
{
const
double
collisionProbability
=
1.0
f
-
exp
(
-
collisionFrequency
*
stepSize
);
for
(
int
i
=
0
;
i
<
(
int
)
atomGroups
.
size
();
++
i
)
{
for
(
auto
&
group
:
atomGroups
)
{
if
(
SimTKOpenMMUtilities
::
getUniformlyDistributedRandomNumber
()
<
collisionProbability
)
{
// A collision occurred, so set the velocities to new values chosen from a Boltzmann distribution.
for
(
int
j
=
0
;
j
<
(
int
)
atomGroups
[
i
].
size
();
j
++
)
{
int
atom
=
atomGroups
[
i
][
j
];
for
(
int
atom
:
group
)
{
if
(
atomMasses
[
atom
]
!=
0
)
{
const
double
velocityScale
=
static_cast
<
double
>
(
sqrt
(
BOLTZ
*
temperature
/
atomMasses
[
atom
]));
atomVelocities
[
atom
][
0
]
=
velocityScale
*
SimTKOpenMMUtilities
::
getNormallyDistributedRandomNumber
();
...
...
platforms/reference/src/SimTKReference/ReferenceCCMAAlgorithm.cpp
View file @
28b79b2f
...
...
@@ -145,8 +145,7 @@ ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm(int numberOfAtoms,
vector
<
double
>
matrixValue
;
for
(
int
i
=
0
;
i
<
numberOfConstraints
;
i
++
)
{
matrixRowStart
.
push_back
(
matrixValue
.
size
());
for
(
int
j
=
0
;
j
<
(
int
)
matrix
[
i
].
size
();
j
++
)
{
pair
<
int
,
double
>
element
=
matrix
[
i
][
j
];
for
(
auto
&
element
:
matrix
[
i
])
{
matrixColIndex
.
push_back
(
element
.
first
);
matrixValue
.
push_back
(
element
.
second
);
}
...
...
@@ -292,10 +291,8 @@ void ReferenceCCMAAlgorithm::applyConstraints(vector<Vec3>& atomCoordinates,
if
(
_matrix
.
size
()
>
0
)
{
for
(
int
i
=
0
;
i
<
_numberOfConstraints
;
i
++
)
{
double
sum
=
0.0
;
for
(
int
j
=
0
;
j
<
(
int
)
_matrix
[
i
].
size
();
j
++
)
{
pair
<
int
,
double
>
element
=
_matrix
[
i
][
j
];
for
(
auto
&
element
:
_matrix
[
i
])
sum
+=
element
.
second
*
constraintDelta
[
element
.
first
];
}
tempDelta
[
i
]
=
sum
;
}
constraintDelta
=
tempDelta
;
...
...
platforms/reference/src/SimTKReference/ReferenceCustomAngleIxn.cpp
View file @
28b79b2f
...
...
@@ -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 @
28b79b2f
...
...
@@ -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 @
28b79b2f
...
...
@@ -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 @
28b79b2f
...
...
@@ -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 @
28b79b2f
...
...
@@ -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 @
28b79b2f
...
...
@@ -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 @
28b79b2f
...
...
@@ -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 @
28b79b2f
...
...
@@ -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 @
28b79b2f
...
...
@@ -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 @
28b79b2f
...
...
@@ -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
);
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 @
28b79b2f
...
...
@@ -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 @
28b79b2f
...
...
@@ -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
;
...
...
@@ -581,11 +580,9 @@ 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
++
)
{
// loop over atom pairs
...
...
platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp
View file @
28b79b2f
...
...
@@ -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
;
}
}
...
...
Prev
1
2
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