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
083bc501
"vscode:/vscode.git/clone" did not exist on "bf3a4db65fc57f8bbb93726ae8e649922aa98372"
Commit
083bc501
authored
Mar 24, 2017
by
peastman
Browse files
Use C++11 style loops
parent
88313407
Changes
18
Hide whitespace changes
Inline
Side-by-side
Showing
18 changed files
with
172 additions
and
202 deletions
+172
-202
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
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp
.../reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp
+4
-6
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaExtrapolatedPolarization.cpp
...nce/tests/TestReferenceAmoebaExtrapolatedPolarization.cpp
+2
-2
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaGeneralizedKirkwoodForce.cpp
...nce/tests/TestReferenceAmoebaGeneralizedKirkwoodForce.cpp
+2
-2
plugins/amoeba/serialization/src/AmoebaBondForceProxy.cpp
plugins/amoeba/serialization/src/AmoebaBondForceProxy.cpp
+1
-3
plugins/amoeba/serialization/src/AmoebaStretchBendForceProxy.cpp
.../amoeba/serialization/src/AmoebaStretchBendForceProxy.cpp
+1
-2
plugins/cpupme/src/CpuPmeKernels.cpp
plugins/cpupme/src/CpuPmeKernels.cpp
+8
-8
plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp
...s/drude/platforms/reference/src/ReferenceDrudeKernels.cpp
+1
-2
plugins/drude/serialization/src/DrudeForceProxy.cpp
plugins/drude/serialization/src/DrudeForceProxy.cpp
+2
-6
plugins/rpmd/openmmapi/src/RPMDIntegrator.cpp
plugins/rpmd/openmmapi/src/RPMDIntegrator.cpp
+8
-9
plugins/rpmd/platforms/cuda/src/CudaRpmdKernels.cpp
plugins/rpmd/platforms/cuda/src/CudaRpmdKernels.cpp
+8
-8
plugins/rpmd/platforms/opencl/src/OpenCLRpmdKernels.cpp
plugins/rpmd/platforms/opencl/src/OpenCLRpmdKernels.cpp
+10
-10
plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.cpp
...ins/rpmd/platforms/reference/src/ReferenceRpmdKernels.cpp
+9
-9
No files found.
plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp
View file @
083bc501
...
@@ -161,14 +161,14 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
...
@@ -161,14 +161,14 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
// Double loop over different atom types.
// Double loop over different atom types.
std
::
string
sigmaCombiningRule
=
force
.
getSigmaCombiningRule
();
std
::
string
sigmaCombiningRule
=
force
.
getSigmaCombiningRule
();
std
::
string
epsilonCombiningRule
=
force
.
getEpsilonCombiningRule
();
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
;
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.
// AMOEBA combining rules, copied over from the CUDA code.
double
iSigma
=
class1
->
first
.
first
;
double
iSigma
=
class1
.
first
.
first
;
double
jSigma
=
class2
->
first
.
first
;
double
jSigma
=
class2
.
first
.
first
;
double
iEpsilon
=
class1
->
first
.
second
;
double
iEpsilon
=
class1
.
first
.
second
;
double
jEpsilon
=
class2
->
first
.
second
;
double
jEpsilon
=
class2
.
first
.
second
;
// ARITHMETIC = 1
// ARITHMETIC = 1
// GEOMETRIC = 2
// GEOMETRIC = 2
// CUBIC-MEAN = 3
// CUBIC-MEAN = 3
...
@@ -207,7 +207,7 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
...
@@ -207,7 +207,7 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
epsilon
=
0.0
;
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.
// Below is an exact copy of stuff from the previous block.
double
rv
=
sigma
;
double
rv
=
sigma
;
double
termik
=
2.0
*
M_PI
*
count
;
// termik is equivalent to 2 * pi * count.
double
termik
=
2.0
*
M_PI
*
count
;
// termik is equivalent to 2 * pi * count.
...
...
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @
083bc501
...
@@ -973,8 +973,8 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
...
@@ -973,8 +973,8 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
molecularQuadrupolesVec
.
push_back
((
float
)
quadrupole
[
5
]);
molecularQuadrupolesVec
.
push_back
((
float
)
quadrupole
[
5
]);
}
}
hasQuadrupoles
=
false
;
hasQuadrupoles
=
false
;
for
(
int
i
=
0
;
i
<
(
int
)
molecularQuadrupolesVec
.
size
();
i
++
)
for
(
auto
q
:
molecularQuadrupolesVec
)
if
(
molecularQuadrupolesVec
[
i
]
!=
0.0
)
if
(
q
!=
0.0
)
hasQuadrupoles
=
true
;
hasQuadrupoles
=
true
;
int
paddedNumAtoms
=
cu
.
getPaddedNumAtoms
();
int
paddedNumAtoms
=
cu
.
getPaddedNumAtoms
();
for
(
int
i
=
numMultipoles
;
i
<
paddedNumAtoms
;
i
++
)
{
for
(
int
i
=
numMultipoles
;
i
<
paddedNumAtoms
;
i
++
)
{
...
@@ -1049,15 +1049,15 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
...
@@ -1049,15 +1049,15 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
Covalent13
,
atoms
);
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
Covalent13
,
atoms
);
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
for
(
set
<
int
>::
const_iterator
iter
=
allAtoms
.
begin
();
iter
!=
allAtoms
.
end
();
++
iter
)
for
(
int
atom
:
allAtoms
)
covalentFlagValues
.
push_back
(
make_int3
(
i
,
*
iter
,
0
));
covalentFlagValues
.
push_back
(
make_int3
(
i
,
atom
,
0
));
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
Covalent14
,
atoms
);
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
Covalent14
,
atoms
);
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
for
(
int
j
=
0
;
j
<
(
int
)
atoms
.
size
();
j
++
)
for
(
int
atom
:
atoms
)
covalentFlagValues
.
push_back
(
make_int3
(
i
,
atom
s
[
j
]
,
1
));
covalentFlagValues
.
push_back
(
make_int3
(
i
,
atom
,
1
));
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
Covalent15
,
atoms
);
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
Covalent15
,
atoms
);
for
(
int
j
=
0
;
j
<
(
int
)
atoms
.
size
();
j
++
)
for
(
int
atom
:
atoms
)
covalentFlagValues
.
push_back
(
make_int3
(
i
,
atom
s
[
j
]
,
2
));
covalentFlagValues
.
push_back
(
make_int3
(
i
,
atom
,
2
));
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
PolarizationCovalent11
,
atoms
);
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
PolarizationCovalent11
,
atoms
);
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
allAtoms
.
insert
(
atoms
.
begin
(),
atoms
.
end
());
...
@@ -1068,15 +1068,14 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
...
@@ -1068,15 +1068,14 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
vector
<
int
>
atoms12
;
vector
<
int
>
atoms12
;
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
PolarizationCovalent12
,
atoms12
);
force
.
getCovalentMap
(
i
,
AmoebaMultipoleForce
::
PolarizationCovalent12
,
atoms12
);
for
(
int
j
=
0
;
j
<
(
int
)
atoms
.
size
();
j
++
)
for
(
int
atom
:
atoms
)
if
(
find
(
atoms12
.
begin
(),
atoms12
.
end
(),
atom
s
[
j
]
)
==
atoms12
.
end
())
if
(
find
(
atoms12
.
begin
(),
atoms12
.
end
(),
atom
)
==
atoms12
.
end
())
polarizationFlagValues
.
push_back
(
make_int2
(
i
,
atom
s
[
j
]
));
polarizationFlagValues
.
push_back
(
make_int2
(
i
,
atom
));
}
}
set
<
pair
<
int
,
int
>
>
tilesWithExclusions
;
set
<
pair
<
int
,
int
>
>
tilesWithExclusions
;
for
(
int
atom1
=
0
;
atom1
<
(
int
)
exclusions
.
size
();
++
atom1
)
{
for
(
int
atom1
=
0
;
atom1
<
(
int
)
exclusions
.
size
();
++
atom1
)
{
int
x
=
atom1
/
CudaContext
::
TileSize
;
int
x
=
atom1
/
CudaContext
::
TileSize
;
for
(
int
j
=
0
;
j
<
(
int
)
exclusions
[
atom1
].
size
();
++
j
)
{
for
(
int
atom2
:
exclusions
[
atom1
])
{
int
atom2
=
exclusions
[
atom1
][
j
];
int
y
=
atom2
/
CudaContext
::
TileSize
;
int
y
=
atom2
/
CudaContext
::
TileSize
;
tilesWithExclusions
.
insert
(
make_pair
(
max
(
x
,
y
),
min
(
x
,
y
)));
tilesWithExclusions
.
insert
(
make_pair
(
max
(
x
,
y
),
min
(
x
,
y
)));
}
}
...
@@ -1412,10 +1411,10 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
...
@@ -1412,10 +1411,10 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
}
}
covalentFlags
=
CudaArray
::
create
<
uint2
>
(
cu
,
nb
.
getExclusions
().
getSize
(),
"covalentFlags"
);
covalentFlags
=
CudaArray
::
create
<
uint2
>
(
cu
,
nb
.
getExclusions
().
getSize
(),
"covalentFlags"
);
vector
<
uint2
>
covalentFlagsVec
(
nb
.
getExclusions
().
getSize
(),
make_uint2
(
0
,
0
));
vector
<
uint2
>
covalentFlagsVec
(
nb
.
getExclusions
().
getSize
(),
make_uint2
(
0
,
0
));
for
(
int
i
=
0
;
i
<
(
int
)
covalentFlagValues
.
size
();
i
++
)
{
for
(
int
3
values
:
covalentFlagValues
)
{
int
atom1
=
co
val
entFlagValues
[
i
]
.
x
;
int
atom1
=
val
ues
.
x
;
int
atom2
=
co
val
entFlagValues
[
i
]
.
y
;
int
atom2
=
val
ues
.
y
;
int
value
=
co
val
entFlagValues
[
i
]
.
z
;
int
value
=
val
ues
.
z
;
int
x
=
atom1
/
CudaContext
::
TileSize
;
int
x
=
atom1
/
CudaContext
::
TileSize
;
int
offset1
=
atom1
-
x
*
CudaContext
::
TileSize
;
int
offset1
=
atom1
-
x
*
CudaContext
::
TileSize
;
int
y
=
atom2
/
CudaContext
::
TileSize
;
int
y
=
atom2
/
CudaContext
::
TileSize
;
...
@@ -1446,9 +1445,9 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
...
@@ -1446,9 +1445,9 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
polarizationGroupFlags
=
CudaArray
::
create
<
unsigned
int
>
(
cu
,
nb
.
getExclusions
().
getSize
(),
"polarizationGroupFlags"
);
polarizationGroupFlags
=
CudaArray
::
create
<
unsigned
int
>
(
cu
,
nb
.
getExclusions
().
getSize
(),
"polarizationGroupFlags"
);
vector
<
unsigned
int
>
polarizationGroupFlagsVec
(
nb
.
getExclusions
().
getSize
(),
0
);
vector
<
unsigned
int
>
polarizationGroupFlagsVec
(
nb
.
getExclusions
().
getSize
(),
0
);
for
(
int
i
=
0
;
i
<
(
int
)
polarizationFlagValues
.
size
();
i
++
)
{
for
(
int
2
values
:
polarizationFlagValues
)
{
int
atom1
=
polarizationFlagV
alues
[
i
]
.
x
;
int
atom1
=
v
alues
.
x
;
int
atom2
=
polarizationFlagV
alues
[
i
]
.
y
;
int
atom2
=
v
alues
.
y
;
int
x
=
atom1
/
CudaContext
::
TileSize
;
int
x
=
atom1
/
CudaContext
::
TileSize
;
int
offset1
=
atom1
-
x
*
CudaContext
::
TileSize
;
int
offset1
=
atom1
-
x
*
CudaContext
::
TileSize
;
int
y
=
atom2
/
CudaContext
::
TileSize
;
int
y
=
atom2
/
CudaContext
::
TileSize
;
...
@@ -1473,10 +1472,12 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
...
@@ -1473,10 +1472,12 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
double
CudaCalcAmoebaMultipoleForceKernel
::
execute
(
ContextImpl
&
context
,
bool
includeForces
,
bool
includeEnergy
)
{
double
CudaCalcAmoebaMultipoleForceKernel
::
execute
(
ContextImpl
&
context
,
bool
includeForces
,
bool
includeEnergy
)
{
if
(
!
hasInitializedScaleFactors
)
{
if
(
!
hasInitializedScaleFactors
)
{
initializeScaleFactors
();
initializeScaleFactors
();
for
(
int
i
=
0
;
i
<
(
int
)
context
.
getForceImpls
()
.
size
()
&&
gkKernel
==
NULL
;
i
++
)
{
for
(
auto
impl
:
context
.
getForceImpls
())
{
AmoebaGeneralizedKirkwoodForceImpl
*
gkImpl
=
dynamic_cast
<
AmoebaGeneralizedKirkwoodForceImpl
*>
(
context
.
getForceImpls
()[
i
]
);
AmoebaGeneralizedKirkwoodForceImpl
*
gkImpl
=
dynamic_cast
<
AmoebaGeneralizedKirkwoodForceImpl
*>
(
impl
);
if
(
gkImpl
!=
NULL
)
if
(
gkImpl
!=
NULL
)
{
gkKernel
=
dynamic_cast
<
CudaCalcAmoebaGeneralizedKirkwoodForceKernel
*>
(
&
gkImpl
->
getKernel
().
getImpl
());
gkKernel
=
dynamic_cast
<
CudaCalcAmoebaGeneralizedKirkwoodForceKernel
*>
(
&
gkImpl
->
getKernel
().
getImpl
());
break
;
}
}
}
}
}
CudaNonbondedUtilities
&
nb
=
cu
.
getNonbondedUtilities
();
CudaNonbondedUtilities
&
nb
=
cu
.
getNonbondedUtilities
();
...
@@ -2232,8 +2233,8 @@ void CudaCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImpl& co
...
@@ -2232,8 +2233,8 @@ void CudaCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImpl& co
molecularQuadrupolesVec
.
push_back
((
float
)
quadrupole
[
5
]);
molecularQuadrupolesVec
.
push_back
((
float
)
quadrupole
[
5
]);
}
}
if
(
!
hasQuadrupoles
)
{
if
(
!
hasQuadrupoles
)
{
for
(
int
i
=
0
;
i
<
(
int
)
molecularQuadrupolesVec
.
size
();
i
++
)
for
(
auto
q
:
molecularQuadrupolesVec
)
if
(
molecularQuadrupolesVec
[
i
]
!=
0.0
)
if
(
q
!=
0.0
)
throw
OpenMMException
(
"updateParametersInContext: Cannot set a non-zero quadrupole moment, because quadrupoles were excluded from the kernel"
);
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
++
)
{
for
(
int
i
=
force
.
getNumMultipoles
();
i
<
cu
.
getPaddedNumAtoms
();
i
++
)
{
...
...
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaExtrapolatedPolarization.cpp
View file @
083bc501
...
@@ -290,8 +290,8 @@ static void check_finite_differences(vector<Vec3> analytic_forces, Context &cont
...
@@ -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.
// 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
;
double
norm
=
0.0
;
for
(
int
i
=
0
;
i
<
(
int
)
analytic_forces
.
size
();
++
i
)
for
(
auto
&
f
:
analytic_forces
)
norm
+=
analytic_forces
[
i
].
dot
(
analytic_forces
[
i
]
);
norm
+=
f
.
dot
(
f
);
norm
=
std
::
sqrt
(
norm
);
norm
=
std
::
sqrt
(
norm
);
const
double
stepSize
=
1e-3
;
const
double
stepSize
=
1e-3
;
double
step
=
0.5
*
stepSize
/
norm
;
double
step
=
0.5
*
stepSize
/
norm
;
...
...
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp
View file @
083bc501
...
@@ -61,8 +61,8 @@ static void checkFiniteDifferences(vector<Vec3> forces, Context &context, vector
...
@@ -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.
// 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;
double norm = 0.0;
for (
int i = 0; i < (int) forces.size(); ++i
)
for (
auto& f : forces
)
norm += f
orces[i].dot(forces[i]
);
norm += f
.dot(f
);
norm = std::sqrt(norm);
norm = std::sqrt(norm);
const double stepSize = 1e-3;
const double stepSize = 1e-3;
double step = 0.5*stepSize/norm;
double step = 0.5*stepSize/norm;
...
...
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaMultipoleForce.cpp
View file @
083bc501
...
@@ -2919,8 +2919,8 @@ static void testNoQuadrupoles(bool usePme) {
...
@@ -2919,8 +2919,8 @@ static void testNoQuadrupoles(bool usePme) {
int
axisType
,
atomX
,
atomY
,
atomZ
;
int
axisType
,
atomX
,
atomY
,
atomZ
;
vector
<
double
>
dipole
,
quadrupole
;
vector
<
double
>
dipole
,
quadrupole
;
amoebaMultipoleForce
->
getMultipoleParameters
(
i
,
charge
,
dipole
,
quadrupole
,
axisType
,
atomZ
,
atomX
,
atomY
,
thole
,
damping
,
polarity
);
amoebaMultipoleForce
->
getMultipoleParameters
(
i
,
charge
,
dipole
,
quadrupole
,
axisType
,
atomZ
,
atomX
,
atomY
,
thole
,
damping
,
polarity
);
for
(
int
j
=
0
;
j
<
(
int
)
quadrupole
.
size
();
j
++
)
for
(
auto
&
q
:
quadrupole
)
q
uadrupole
[
j
]
=
0
;
q
=
0
;
amoebaMultipoleForce
->
setMultipoleParameters
(
i
,
charge
,
dipole
,
quadrupole
,
axisType
,
atomZ
,
atomX
,
atomY
,
thole
,
damping
,
polarity
);
amoebaMultipoleForce
->
setMultipoleParameters
(
i
,
charge
,
dipole
,
quadrupole
,
axisType
,
atomZ
,
atomX
,
atomY
,
thole
,
damping
,
polarity
);
}
}
amoebaMultipoleForce
->
updateParametersInContext
(
context
);
amoebaMultipoleForce
->
updateParametersInContext
(
context
);
...
...
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp
View file @
083bc501
...
@@ -206,18 +206,18 @@ void AmoebaReferenceMultipoleForce::setupScaleMaps(const vector< vector< vector<
...
@@ -206,18 +206,18 @@ void AmoebaReferenceMultipoleForce::setupScaleMaps(const vector< vector< vector<
// pScale & mScale
// pScale & mScale
for (unsigned jj = 0; jj < AmoebaMultipoleForce::PolarizationCovalent11; jj++) {
for (unsigned jj = 0; jj < AmoebaMultipoleForce::PolarizationCovalent11; jj++) {
const
vector
<
int
>
covalentList
=
covalentInfo
[
jj
];
for (int covalentIndex : covalentInfo[jj]) {
for
(
unsigned
int
kk
=
0
;
kk
<
covalentList
.
size
();
kk
++
)
{
if (covalentIndex < ii)
unsigned
int
covalentIndex
=
static_cast
<
unsigned
int
>
(
covalentList
[
kk
]);
continue;
if
(
covalentIndex
<
ii
)
continue
;
// handle 0.5 factor for p14
// handle 0.5 factor for p14
int hit = 0;
int hit = 0;
if (jj == AmoebaMultipoleForce::Covalent14) {
if (jj == AmoebaMultipoleForce::Covalent14) {
for
(
unsigned
int
mm
=
0
;
mm
<
covalentListP11
.
size
()
&&
hit
==
0
;
mm
++
)
{
for (
int mm : covalentListP11
) {
if
(
covalentListP11
[
mm
]
==
covalentIndex
)
{
if (
mm
== covalentIndex) {
hit = 1;
hit = 1;
break;
}
}
}
}
}
}
...
@@ -231,9 +231,7 @@ void AmoebaReferenceMultipoleForce::setupScaleMaps(const vector< vector< vector<
...
@@ -231,9 +231,7 @@ void AmoebaReferenceMultipoleForce::setupScaleMaps(const vector< vector< vector<
// dScale & uScale
// dScale & uScale
for (unsigned jj = AmoebaMultipoleForce::PolarizationCovalent11; jj < covalentInfo.size(); jj++) {
for (unsigned jj = AmoebaMultipoleForce::PolarizationCovalent11; jj < covalentInfo.size(); jj++) {
const
vector
<
int
>
covalentList
=
covalentInfo
[
jj
];
for (int covalentIndex : covalentInfo[jj]) {
for
(
unsigned
int
kk
=
0
;
kk
<
covalentList
.
size
();
kk
++
)
{
unsigned
int
covalentIndex
=
static_cast
<
unsigned
int
>
(
covalentList
[
kk
]);
if (covalentIndex < ii)continue;
if (covalentIndex < ii)continue;
_scaleMaps[ii][D_SCALE][covalentIndex] = _dScale[jj-4];
_scaleMaps[ii][D_SCALE][covalentIndex] = _dScale[jj-4];
_scaleMaps[ii][U_SCALE][covalentIndex] = _uScale[jj-4];
_scaleMaps[ii][U_SCALE][covalentIndex] = _uScale[jj-4];
...
@@ -827,16 +825,16 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
...
@@ -827,16 +825,16 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
double rr3 = -rrI[0];
double rr3 = -rrI[0];
double rr5 = rrI[1];
double rr5 = rrI[1];
for
(
unsigned
int
ii
=
0
;
ii
<
updateInducedDipoleFields
.
size
();
ii
++
)
{
for (
auto& field :
updateInducedDipoleFields) {
calculateInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, rr3, rr5, deltaR,
calculateInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, rr3, rr5, deltaR,
*
updateInducedDipoleFields
[
ii
].
inducedDipoles
,
updateInducedDipoleFields
[
ii
]
.
inducedDipoleField
);
*
field.inducedDipoles, field
.inducedDipoleField);
if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated) {
if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated) {
// Compute and store the field gradient for later use.
// Compute and store the field gradient for later use.
double dx = deltaR[0];
double dx = deltaR[0];
double dy = deltaR[1];
double dy = deltaR[1];
double dz = deltaR[2];
double dz = deltaR[2];
OpenMM
::
Vec3
&
dipolesI
=
(
*
updateInducedDipoleFields
[
ii
]
.
inducedDipoles
)[
particleI
.
particleIndex
];
OpenMM::Vec3 &dipolesI = (*
field
.inducedDipoles)[particleI.particleIndex];
double xDipole = dipolesI[0];
double xDipole = dipolesI[0];
double yDipole = dipolesI[1];
double yDipole = dipolesI[1];
double zDipole = dipolesI[2];
double zDipole = dipolesI[2];
...
@@ -848,14 +846,14 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
...
@@ -848,14 +846,14 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
double Exz = muDotR*dx*dz*rrI[2] - (xDipole*dz + zDipole*dx)*rrI[1];
double Exz = muDotR*dx*dz*rrI[2] - (xDipole*dz + zDipole*dx)*rrI[1];
double Eyz = muDotR*dy*dz*rrI[2] - (yDipole*dz + zDipole*dy)*rrI[1];
double Eyz = muDotR*dy*dz*rrI[2] - (yDipole*dz + zDipole*dy)*rrI[1];
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
0
]
-=
Exx
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][0] -= Exx;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
1
]
-=
Eyy
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][1] -= Eyy;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
2
]
-=
Ezz
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][2] -= Ezz;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
3
]
-=
Exy
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][3] -= Exy;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
4
]
-=
Exz
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][4] -= Exz;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
5
]
-=
Eyz
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][5] -= Eyz;
OpenMM
::
Vec3
&
dipolesJ
=
(
*
updateInducedDipoleFields
[
ii
]
.
inducedDipoles
)[
particleJ
.
particleIndex
];
OpenMM::Vec3 &dipolesJ = (*
field
.inducedDipoles)[particleJ.particleIndex];
xDipole = dipolesJ[0];
xDipole = dipolesJ[0];
yDipole = dipolesJ[1];
yDipole = dipolesJ[1];
zDipole = dipolesJ[2];
zDipole = dipolesJ[2];
...
@@ -867,12 +865,12 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
...
@@ -867,12 +865,12 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
Exz = muDotR*dx*dz*rrI[2] - (xDipole*dz + zDipole*dx)*rrI[1];
Exz = muDotR*dx*dz*rrI[2] - (xDipole*dz + zDipole*dx)*rrI[1];
Eyz = muDotR*dy*dz*rrI[2] - (yDipole*dz + zDipole*dy)*rrI[1];
Eyz = muDotR*dy*dz*rrI[2] - (yDipole*dz + zDipole*dy)*rrI[1];
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
0
]
+=
Exx
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][0] += Exx;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
1
]
+=
Eyy
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][1] += Eyy;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
2
]
+=
Ezz
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][2] += Ezz;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
3
]
+=
Exy
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][3] += Exy;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
4
]
+=
Exz
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][4] += Exz;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
5
]
+=
Eyz
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][5] += Eyz;
}
}
}
}
}
}
...
@@ -881,8 +879,8 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipoleFields(const vector<Mu
...
@@ -881,8 +879,8 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipoleFields(const vector<Mu
// Initialize the fields to zero.
// Initialize the fields to zero.
Vec3 zeroVec(0.0, 0.0, 0.0);
Vec3 zeroVec(0.0, 0.0, 0.0);
for
(
unsigned
int
ii
=
0
;
ii
<
updateInducedDipoleFields
.
size
();
ii
++
)
for (
auto& field :
updateInducedDipoleFields)
std
::
fill
(
updateInducedDipoleFields
[
ii
].
inducedDipoleField
.
begin
(),
updateInducedDipoleFields
[
ii
]
.
inducedDipoleField
.
end
(),
zeroVec
);
std::fill(
field.inducedDipoleField.begin(), field
.inducedDipoleField.end(), zeroVec);
// Add fields from all induced dipoles.
// Add fields from all induced dipoles.
...
@@ -901,11 +899,11 @@ double AmoebaReferenceMultipoleForce::updateInducedDipoleFields(const vector<Mul
...
@@ -901,11 +899,11 @@ double AmoebaReferenceMultipoleForce::updateInducedDipoleFields(const vector<Mul
// Update the induced dipoles and calculate the convergence factor, maxEpsilon
// Update the induced dipoles and calculate the convergence factor, maxEpsilon
double maxEpsilon = 0.0;
double maxEpsilon = 0.0;
for
(
unsigned
int
kk
=
0
;
kk
<
updateInducedDipoleFields
.
size
();
kk
++
)
{
for (
auto& field :
updateInducedDipoleFields) {
double epsilon = updateInducedDipole(particleData,
double epsilon = updateInducedDipole(particleData,
*
updateInducedDipoleFields
[
kk
]
.
fixedMultipoleField
,
*
field
.fixedMultipoleField,
updateInducedDipoleFields
[
kk
]
.
inducedDipoleField
,
field
.inducedDipoleField,
*
updateInducedDipoleFields
[
kk
]
.
inducedDipoles
);
*
field
.inducedDipoles);
maxEpsilon = epsilon > maxEpsilon ? epsilon : maxEpsilon;
maxEpsilon = epsilon > maxEpsilon ? epsilon : maxEpsilon;
}
}
...
@@ -1792,9 +1790,8 @@ double AmoebaReferenceMultipoleForce::calculateElectrostatic(const vector<Multip
...
@@ -1792,9 +1790,8 @@ double AmoebaReferenceMultipoleForce::calculateElectrostatic(const vector<Multip
{
{
double energy = 0.0;
double energy = 0.0;
vector<double> scaleFactors(LAST_SCALE_TYPE_INDEX);
vector<double> scaleFactors(LAST_SCALE_TYPE_INDEX);
for
(
unsigned
int
kk
=
0
;
kk
<
scaleFactors
.
size
();
kk
++
)
{
for (auto& s : scaleFactors)
scaleFactors
[
kk
]
=
1.0
;
s = 1.0;
}
// main loop over particle pairs
// main loop over particle pairs
...
@@ -1971,9 +1968,7 @@ void AmoebaReferenceMultipoleForce::calculateLabFramePermanentDipoles(const vect
...
@@ -1971,9 +1968,7 @@ void AmoebaReferenceMultipoleForce::calculateLabFramePermanentDipoles(const vect
multipoleAtomCovalentInfo, particleData);
multipoleAtomCovalentInfo, particleData);
outputRotatedPermanentDipoles.resize(_numParticles);
outputRotatedPermanentDipoles.resize(_numParticles);
for (int i = 0; i < _numParticles; i++)
for (int i = 0; i < _numParticles; i++)
{
outputRotatedPermanentDipoles[i] = particleData[i].dipole;
outputRotatedPermanentDipoles
[
i
]
=
particleData
[
i
].
dipole
;
}
}
}
void AmoebaReferenceMultipoleForce::calculateTotalDipoles(const vector<Vec3>& particlePositions,
void AmoebaReferenceMultipoleForce::calculateTotalDipoles(const vector<Vec3>& particlePositions,
...
@@ -1997,12 +1992,8 @@ void AmoebaReferenceMultipoleForce::calculateTotalDipoles(const vector<Vec3>& pa
...
@@ -1997,12 +1992,8 @@ void AmoebaReferenceMultipoleForce::calculateTotalDipoles(const vector<Vec3>& pa
multipoleAtomCovalentInfo, particleData);
multipoleAtomCovalentInfo, particleData);
outputTotalDipoles.resize(_numParticles);
outputTotalDipoles.resize(_numParticles);
for (int i = 0; i < _numParticles; i++)
for (int i = 0; i < _numParticles; i++)
{
for (int j = 0; j < 3; j++)
for
(
int
j
=
0
;
j
<
3
;
j
++
)
outputTotalDipoles[i][j] = particleData[i].dipole[j] + _inducedDipole[i][j];
{
outputTotalDipoles
[
i
][
j
]
=
particleData
[
i
].
dipole
[
j
]
+
_inducedDipole
[
i
][
j
];
}
}
}
}
void AmoebaReferenceMultipoleForce::calculateAmoebaSystemMultipoleMoments(const vector<double>& masses,
void AmoebaReferenceMultipoleForce::calculateAmoebaSystemMultipoleMoments(const vector<double>& masses,
...
@@ -2173,9 +2164,8 @@ void AmoebaReferenceMultipoleForce::calculateElectrostaticPotential(const vector
...
@@ -2173,9 +2164,8 @@ void AmoebaReferenceMultipoleForce::calculateElectrostaticPotential(const vector
multipoleAtomCovalentInfo, particleData);
multipoleAtomCovalentInfo, particleData);
potential.resize(grid.size());
potential.resize(grid.size());
for
(
unsigned
int
ii
=
0
;
ii
<
grid
.
size
();
ii
++
)
{
for (auto& p : potential)
potential
[
ii
]
=
0.0
;
p = 0.0;
}
for (unsigned int ii = 0; ii < _numParticles; ii++) {
for (unsigned int ii = 0; ii < _numParticles; ii++) {
for (unsigned int jj = 0; jj < grid.size(); jj++) {
for (unsigned int jj = 0; jj < grid.size(); jj++) {
...
@@ -2184,9 +2174,8 @@ void AmoebaReferenceMultipoleForce::calculateElectrostaticPotential(const vector
...
@@ -2184,9 +2174,8 @@ void AmoebaReferenceMultipoleForce::calculateElectrostaticPotential(const vector
}
}
double term = _electric/_dielectric;
double term = _electric/_dielectric;
for
(
unsigned
int
ii
=
0
;
ii
<
grid
.
size
();
ii
++
)
{
for (auto& p : potential)
potential
[
ii
]
*=
term
;
p *= term;
}
}
}
AmoebaReferenceMultipoleForce::UpdateInducedDipoleFieldStruct::UpdateInducedDipoleFieldStruct(vector<OpenMM::Vec3>& inputFixed_E_Field, vector<OpenMM::Vec3>& inputInducedDipoles, vector<vector<Vec3> >& extrapolatedDipoles, vector<vector<double> >& extrapolatedDipoleFieldGradient) :
AmoebaReferenceMultipoleForce::UpdateInducedDipoleFieldStruct::UpdateInducedDipoleFieldStruct(vector<OpenMM::Vec3>& inputFixed_E_Field, vector<OpenMM::Vec3>& inputInducedDipoles, vector<vector<Vec3> >& extrapolatedDipoles, vector<vector<double> >& extrapolatedDipoleFieldGradient) :
...
@@ -4066,9 +4055,8 @@ double AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateElectrostatic(
...
@@ -4066,9 +4055,8 @@ double AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateElectrostatic(
// correct vacuum to SCRF derivatives (ediff1 in TINKER)
// correct vacuum to SCRF derivatives (ediff1 in TINKER)
vector<double> scaleFactors(LAST_SCALE_TYPE_INDEX);
vector<double> scaleFactors(LAST_SCALE_TYPE_INDEX);
for
(
unsigned
int
kk
=
0
;
kk
<
scaleFactors
.
size
();
kk
++
)
{
for (auto& s : scaleFactors)
scaleFactors
[
kk
]
=
1.0
;
s = 1.0;
}
double eDiffEnergy = 0.0;
double eDiffEnergy = 0.0;
for (unsigned int ii = 0; ii < particleData.size(); ii++) {
for (unsigned int ii = 0; ii < particleData.size(); ii++) {
...
@@ -4082,9 +4070,8 @@ double AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateElectrostatic(
...
@@ -4082,9 +4070,8 @@ double AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateElectrostatic(
scaleFactors[P_SCALE], scaleFactors[D_SCALE], forces, torques);
scaleFactors[P_SCALE], scaleFactors[D_SCALE], forces, torques);
if (jj <= _maxScaleIndex[ii]) {
if (jj <= _maxScaleIndex[ii]) {
for
(
unsigned
int
kk
=
0
;
kk
<
LAST_SCALE_TYPE_INDEX
;
kk
++
)
{
for (auto& s : scaleFactors)
scaleFactors
[
kk
]
=
1.0
;
s = 1.0;
}
}
}
}
}
}
}
...
@@ -5012,9 +4999,8 @@ void AmoebaReferencePmeMultipoleForce::initializePmeGrid()
...
@@ -5012,9 +4999,8 @@ void AmoebaReferencePmeMultipoleForce::initializePmeGrid()
if (_pmeGrid == NULL)
if (_pmeGrid == NULL)
return;
return;
for
(
int
jj
=
0
;
jj
<
_totalGridSize
;
jj
++
)
{
for (int jj = 0; jj < _totalGridSize; jj++)
_pmeGrid[jj].re = _pmeGrid[jj].im = 0.0;
_pmeGrid[jj].re = _pmeGrid[jj].im = 0.0;
}
}
}
void AmoebaReferencePmeMultipoleForce::getPeriodicDelta(Vec3& deltaR) const
void AmoebaReferencePmeMultipoleForce::getPeriodicDelta(Vec3& deltaR) const
...
@@ -5137,9 +5123,9 @@ void AmoebaReferencePmeMultipoleForce::initializeBSplineModuli()
...
@@ -5137,9 +5123,9 @@ void AmoebaReferencePmeMultipoleForce::initializeBSplineModuli()
if (i > size/2)
if (i > size/2)
k = k - size;
k = k - size;
double zeta;
double zeta;
if
(
k
==
0
)
{
if (k == 0)
zeta = 1.0;
zeta = 1.0;
}
else
{
else {
double sum1 = 1.0;
double sum1 = 1.0;
double sum2 = 1.0;
double sum2 = 1.0;
factor = M_PI*k/size;
factor = M_PI*k/size;
...
@@ -5270,9 +5256,9 @@ void AmoebaReferencePmeMultipoleForce::calculateFixedMultipoleField(const vector
...
@@ -5270,9 +5256,9 @@ void AmoebaReferencePmeMultipoleForce::calculateFixedMultipoleField(const vector
double term = (4.0/3.0)*(_alphaEwald*_alphaEwald*_alphaEwald)/SQRT_PI;
double term = (4.0/3.0)*(_alphaEwald*_alphaEwald*_alphaEwald)/SQRT_PI;
for (unsigned int jj = 0; jj < _numParticles; jj++) {
for (unsigned int jj = 0; jj < _numParticles; jj++) {
Vec3
selfEnergy
=
particleData
[
jj
].
dipole
*
term
;
Vec3 selfEnergy = particleData[jj].dipole*term;
_fixedMultipoleField
[
jj
]
+=
selfEnergy
;
_fixedMultipoleField[jj] += selfEnergy;
_fixedMultipoleFieldPolar
[
jj
]
=
_fixedMultipoleField
[
jj
];
_fixedMultipoleFieldPolar[jj] = _fixedMultipoleField[jj];
}
}
// include direct space fixed multipole fields
// include direct space fixed multipole fields
...
@@ -5307,9 +5293,8 @@ void AmoebaReferencePmeMultipoleForce::computeBSplinePoint(vector<double4>& thet
...
@@ -5307,9 +5293,8 @@ void AmoebaReferencePmeMultipoleForce::computeBSplinePoint(vector<double4>& thet
int k = i - 1;
int k = i - 1;
double denom = 1.0 / k;
double denom = 1.0 / k;
ARRAY(i,i) = denom * w * ARRAY(k,k);
ARRAY(i,i) = denom * w * ARRAY(k,k);
for
(
int
j
=
1
;
j
<=
i
-
2
;
j
++
)
{
for (int j = 1; j <= i-2; j++)
ARRAY(i,i-j) = denom * ((w+j)*ARRAY(k,i-j-1)+(i-j-w)*ARRAY(k,i-j));
ARRAY(i,i-j) = denom * ((w+j)*ARRAY(k,i-j-1)+(i-j-w)*ARRAY(k,i-j));
}
ARRAY(i,1) = denom * (1.0-w) * ARRAY(k,1);
ARRAY(i,1) = denom * (1.0-w) * ARRAY(k,1);
}
}
...
@@ -5351,9 +5336,8 @@ void AmoebaReferencePmeMultipoleForce::computeBSplinePoint(vector<double4>& thet
...
@@ -5351,9 +5336,8 @@ void AmoebaReferencePmeMultipoleForce::computeBSplinePoint(vector<double4>& thet
// copy coefficients from temporary to permanent storage
// copy coefficients from temporary to permanent storage
for
(
int
i
=
1
;
i
<=
AMOEBA_PME_ORDER
;
i
++
)
{
for (int i = 1; i <= AMOEBA_PME_ORDER; i++)
thetai[i-1] = double4(ARRAY(AMOEBA_PME_ORDER,i), ARRAY(AMOEBA_PME_ORDER-1,i), ARRAY(AMOEBA_PME_ORDER-2,i), ARRAY(AMOEBA_PME_ORDER-3,i));
thetai[i-1] = double4(ARRAY(AMOEBA_PME_ORDER,i), ARRAY(AMOEBA_PME_ORDER-1,i), ARRAY(AMOEBA_PME_ORDER-2,i), ARRAY(AMOEBA_PME_ORDER-3,i));
}
}
}
/**
/**
...
@@ -5377,9 +5361,8 @@ void AmoebaReferencePmeMultipoleForce::computeAmoebaBsplines(const vector<Multip
...
@@ -5377,9 +5361,8 @@ void AmoebaReferencePmeMultipoleForce::computeAmoebaBsplines(const vector<Multip
igrid[jj] += igrid[jj] < 0 ? _pmeGridDimensions[jj] : 0;
igrid[jj] += igrid[jj] < 0 ? _pmeGridDimensions[jj] : 0;
vector<double4> thetaiTemp(AMOEBA_PME_ORDER);
vector<double4> thetaiTemp(AMOEBA_PME_ORDER);
computeBSplinePoint(thetaiTemp, w);
computeBSplinePoint(thetaiTemp, w);
for
(
unsigned
int
kk
=
0
;
kk
<
AMOEBA_PME_ORDER
;
kk
++
)
{
for (unsigned int kk = 0; kk < AMOEBA_PME_ORDER; kk++)
_thetai[jj][ii*AMOEBA_PME_ORDER+kk] = thetaiTemp[kk];
_thetai[jj][ii*AMOEBA_PME_ORDER+kk] = thetaiTemp[kk];
}
}
}
// Record the grid point.
// Record the grid point.
...
@@ -6152,8 +6135,8 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector
...
@@ -6152,8 +6135,8 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector
// Initialize the fields to zero.
// Initialize the fields to zero.
Vec3 zeroVec(0.0, 0.0, 0.0);
Vec3 zeroVec(0.0, 0.0, 0.0);
for
(
unsigned
int
ii
=
0
;
ii
<
updateInducedDipoleFields
.
size
();
ii
++
)
for (
auto& field :
updateInducedDipoleFields)
std
::
fill
(
updateInducedDipoleFields
[
ii
].
inducedDipoleField
.
begin
(),
updateInducedDipoleFields
[
ii
]
.
inducedDipoleField
.
end
(),
zeroVec
);
std::fill(
field.inducedDipoleField.begin(), field
.inducedDipoleField.end(), zeroVec);
// Add fields from direct space interactions.
// Add fields from direct space interactions.
...
@@ -6232,11 +6215,11 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector
...
@@ -6232,11 +6215,11 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector
// self ixn
// self ixn
double term = (4.0/3.0)*(_alphaEwald*_alphaEwald*_alphaEwald)/SQRT_PI;
double term = (4.0/3.0)*(_alphaEwald*_alphaEwald*_alphaEwald)/SQRT_PI;
for
(
unsigned
int
ii
=
0
;
ii
<
updateInducedDipoleFields
.
size
();
ii
++
)
{
for (
auto& field :
updateInducedDipoleFields) {
vector
<
Vec3
>&
inducedDipoles
=
*
updateInducedDipoleFields
[
ii
]
.
inducedDipoles
;
vector<Vec3>& inducedDipoles = *
field
.inducedDipoles;
vector
<
Vec3
>&
field
=
updateI
nducedDipoleField
s
[
ii
]
.
inducedDipoleField
;
vector<Vec3>&
i
nducedDipoleField
= field
.inducedDipoleField;
for (unsigned int jj = 0; jj < particleData.size(); jj++) {
for (unsigned int jj = 0; jj < particleData.size(); jj++) {
f
ield
[
jj
]
+=
inducedDipoles
[
jj
]
*
term
;
inducedDipoleF
ield[jj] += inducedDipoles[jj]*term;
}
}
}
}
}
}
...
@@ -6331,17 +6314,16 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons
...
@@ -6331,17 +6314,16 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons
double preFactor2 = bn2 - rr5;
double preFactor2 = bn2 - rr5;
double preFactor3 = bn3 - rr7;
double preFactor3 = bn3 - rr7;
for
(
unsigned
int
ii
=
0
;
ii
<
updateInducedDipoleFields
.
size
();
ii
++
)
{
for (
auto& field :
updateInducedDipoleFields) {
calculateDirectInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, preFactor1, preFactor2, deltaR,
calculateDirectInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, preFactor1, preFactor2, deltaR,
*
updateInducedDipoleFields
[
ii
].
inducedDipoles
,
*field.inducedDipoles, field.inducedDipoleField);
updateInducedDipoleFields
[
ii
].
inducedDipoleField
);
if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated) {
if (getPolarizationType() == AmoebaReferenceMultipoleForce::Extrapolated) {
// Compute and store the field gradient for later use.
// Compute and store the field gradient for later use.
double dx = deltaR[0];
double dx = deltaR[0];
double dy = deltaR[1];
double dy = deltaR[1];
double dz = deltaR[2];
double dz = deltaR[2];
OpenMM
::
Vec3
&
dipolesI
=
(
*
updateInducedDipoleFields
[
ii
]
.
inducedDipoles
)[
particleI
.
particleIndex
];
OpenMM::Vec3 &dipolesI = (*
field
.inducedDipoles)[particleI.particleIndex];
double xDipole = dipolesI[0];
double xDipole = dipolesI[0];
double yDipole = dipolesI[1];
double yDipole = dipolesI[1];
double zDipole = dipolesI[2];
double zDipole = dipolesI[2];
...
@@ -6353,14 +6335,14 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons
...
@@ -6353,14 +6335,14 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons
double Exz = muDotR*dx*dz*preFactor3 - (xDipole*dz + zDipole*dx)*preFactor2;
double Exz = muDotR*dx*dz*preFactor3 - (xDipole*dz + zDipole*dx)*preFactor2;
double Eyz = muDotR*dy*dz*preFactor3 - (yDipole*dz + zDipole*dy)*preFactor2;
double Eyz = muDotR*dy*dz*preFactor3 - (yDipole*dz + zDipole*dy)*preFactor2;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
0
]
-=
Exx
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][0] -= Exx;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
1
]
-=
Eyy
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][1] -= Eyy;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
2
]
-=
Ezz
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][2] -= Ezz;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
3
]
-=
Exy
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][3] -= Exy;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
4
]
-=
Exz
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][4] -= Exz;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleJ
.
particleIndex
][
5
]
-=
Eyz
;
field
.inducedDipoleFieldGradient[particleJ.particleIndex][5] -= Eyz;
OpenMM
::
Vec3
&
dipolesJ
=
(
*
updateInducedDipoleFields
[
ii
]
.
inducedDipoles
)[
particleJ
.
particleIndex
];
OpenMM::Vec3 &dipolesJ = (*
field
.inducedDipoles)[particleJ.particleIndex];
xDipole = dipolesJ[0];
xDipole = dipolesJ[0];
yDipole = dipolesJ[1];
yDipole = dipolesJ[1];
zDipole = dipolesJ[2];
zDipole = dipolesJ[2];
...
@@ -6372,12 +6354,12 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons
...
@@ -6372,12 +6354,12 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons
Exz = muDotR*dx*dz*preFactor3 - (xDipole*dz + zDipole*dx)*preFactor2;
Exz = muDotR*dx*dz*preFactor3 - (xDipole*dz + zDipole*dx)*preFactor2;
Eyz = muDotR*dy*dz*preFactor3 - (yDipole*dz + zDipole*dy)*preFactor2;
Eyz = muDotR*dy*dz*preFactor3 - (yDipole*dz + zDipole*dy)*preFactor2;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
0
]
+=
Exx
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][0] += Exx;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
1
]
+=
Eyy
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][1] += Eyy;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
2
]
+=
Ezz
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][2] += Ezz;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
3
]
+=
Exy
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][3] += Exy;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
4
]
+=
Exz
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][4] += Exz;
updateInducedDipoleFields
[
ii
]
.
inducedDipoleFieldGradient
[
particleI
.
particleIndex
][
5
]
+=
Eyz
;
field
.inducedDipoleFieldGradient[particleI.particleIndex][5] += Eyz;
}
}
}
}
}
}
...
@@ -6846,9 +6828,8 @@ double AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector<Mul
...
@@ -6846,9 +6828,8 @@ double AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector<Mul
{
{
double energy = 0.0;
double energy = 0.0;
vector<double> scaleFactors(LAST_SCALE_TYPE_INDEX);
vector<double> scaleFactors(LAST_SCALE_TYPE_INDEX);
for
(
unsigned
int
kk
=
0
;
kk
<
scaleFactors
.
size
();
kk
++
)
{
for (auto& s : scaleFactors)
scaleFactors
[
kk
]
=
1.0
;
s = 1.0;
}
// loop over particle pairs for direct space interactions
// loop over particle pairs for direct space interactions
...
@@ -6862,9 +6843,8 @@ double AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector<Mul
...
@@ -6862,9 +6843,8 @@ double AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector<Mul
energy += calculatePmeDirectElectrostaticPairIxn(particleData[ii], particleData[jj], scaleFactors, forces, torques);
energy += calculatePmeDirectElectrostaticPairIxn(particleData[ii], particleData[jj], scaleFactors, forces, torques);
if (jj <= _maxScaleIndex[ii]) {
if (jj <= _maxScaleIndex[ii]) {
for
(
unsigned
int
kk
=
0
;
kk
<
LAST_SCALE_TYPE_INDEX
;
kk
++
)
{
for (auto& s : scaleFactors)
scaleFactors
[
kk
]
=
1.0
;
s = 1.0;
}
}
}
}
}
}
}
...
...
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp
View file @
083bc501
...
@@ -277,9 +277,8 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
...
@@ -277,9 +277,8 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
double
sigmaI
=
sigmas
[
ii
];
double
sigmaI
=
sigmas
[
ii
];
double
epsilonI
=
epsilons
[
ii
];
double
epsilonI
=
epsilons
[
ii
];
for
(
std
::
set
<
int
>::
const_iterator
jj
=
allExclusions
[
ii
].
begin
();
jj
!=
allExclusions
[
ii
].
end
();
jj
++
)
{
for
(
int
jj
:
allExclusions
[
ii
])
exclusions
[
*
jj
]
=
1
;
exclusions
[
jj
]
=
1
;
}
for
(
unsigned
int
jj
=
ii
+
1
;
jj
<
static_cast
<
unsigned
int
>
(
numParticles
);
jj
++
)
{
for
(
unsigned
int
jj
=
ii
+
1
;
jj
<
static_cast
<
unsigned
int
>
(
numParticles
);
jj
++
)
{
if
(
exclusions
[
jj
]
==
0
)
{
if
(
exclusions
[
jj
]
==
0
)
{
...
@@ -310,9 +309,8 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
...
@@ -310,9 +309,8 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
}
}
}
}
for
(
std
::
set
<
int
>::
const_iterator
jj
=
allExclusions
[
ii
].
begin
();
jj
!=
allExclusions
[
ii
].
end
();
jj
++
)
{
for
(
int
jj
:
allExclusions
[
ii
])
exclusions
[
*
jj
]
=
0
;
exclusions
[
jj
]
=
0
;
}
}
}
return
energy
;
return
energy
;
...
...
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaExtrapolatedPolarization.cpp
View file @
083bc501
...
@@ -287,8 +287,8 @@ static void check_finite_differences(vector<Vec3> analytic_forces, Context &cont
...
@@ -287,8 +287,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.
// 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
;
double
norm
=
0.0
;
for
(
int
i
=
0
;
i
<
(
int
)
analytic_forces
.
size
();
++
i
)
for
(
auto
&
f
:
analytic_forces
)
norm
+=
analytic_forces
[
i
].
dot
(
analytic_forces
[
i
]
);
norm
+=
f
.
dot
(
f
);
norm
=
std
::
sqrt
(
norm
);
norm
=
std
::
sqrt
(
norm
);
const
double
stepSize
=
1e-3
;
const
double
stepSize
=
1e-3
;
double
step
=
0.5
*
stepSize
/
norm
;
double
step
=
0.5
*
stepSize
/
norm
;
...
...
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaGeneralizedKirkwoodForce.cpp
View file @
083bc501
...
@@ -61,8 +61,8 @@ static void checkFiniteDifferences(vector<Vec3> forces, Context &context, vector
...
@@ -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.
// 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;
double norm = 0.0;
for (
int i = 0; i < (int) forces.size(); ++i
)
for (
auto& f : forces
)
norm += f
orces[i].dot(forces[i]
);
norm += f
.dot(f
);
norm = std::sqrt(norm);
norm = std::sqrt(norm);
const double stepSize = 1e-3;
const double stepSize = 1e-3;
double step = 0.5*stepSize/norm;
double step = 0.5*stepSize/norm;
...
...
plugins/amoeba/serialization/src/AmoebaBondForceProxy.cpp
View file @
083bc501
...
@@ -72,10 +72,8 @@ void* AmoebaBondForceProxy::deserialize(const SerializationNode& node) const {
...
@@ -72,10 +72,8 @@ void* AmoebaBondForceProxy::deserialize(const SerializationNode& node) const {
force
->
setAmoebaGlobalBondCubic
(
node
.
getDoubleProperty
(
"cubic"
));
force
->
setAmoebaGlobalBondCubic
(
node
.
getDoubleProperty
(
"cubic"
));
force
->
setAmoebaGlobalBondQuartic
(
node
.
getDoubleProperty
(
"quartic"
));
force
->
setAmoebaGlobalBondQuartic
(
node
.
getDoubleProperty
(
"quartic"
));
const
SerializationNode
&
bonds
=
node
.
getChildNode
(
"Bonds"
);
const
SerializationNode
&
bonds
=
node
.
getChildNode
(
"Bonds"
);
for
(
unsigned
int
ii
=
0
;
ii
<
(
int
)
bonds
.
getChildren
().
size
();
ii
++
)
{
for
(
auto
&
bond
:
bonds
.
getChildren
())
const
SerializationNode
&
bond
=
bonds
.
getChildren
()[
ii
];
force
->
addBond
(
bond
.
getIntProperty
(
"p1"
),
bond
.
getIntProperty
(
"p2"
),
bond
.
getDoubleProperty
(
"d"
),
bond
.
getDoubleProperty
(
"k"
));
force
->
addBond
(
bond
.
getIntProperty
(
"p1"
),
bond
.
getIntProperty
(
"p2"
),
bond
.
getDoubleProperty
(
"d"
),
bond
.
getDoubleProperty
(
"k"
));
}
}
}
catch
(...)
{
catch
(...)
{
delete
force
;
delete
force
;
...
...
plugins/amoeba/serialization/src/AmoebaStretchBendForceProxy.cpp
View file @
083bc501
...
@@ -67,8 +67,7 @@ void* AmoebaStretchBendForceProxy::deserialize(const SerializationNode& node) co
...
@@ -67,8 +67,7 @@ void* AmoebaStretchBendForceProxy::deserialize(const SerializationNode& node) co
if
(
version
>
3
)
if
(
version
>
3
)
force
->
setUsesPeriodicBoundaryConditions
(
node
.
getBoolProperty
(
"usesPeriodic"
));
force
->
setUsesPeriodicBoundaryConditions
(
node
.
getBoolProperty
(
"usesPeriodic"
));
const
SerializationNode
&
bonds
=
node
.
getChildNode
(
"StretchBendAngles"
);
const
SerializationNode
&
bonds
=
node
.
getChildNode
(
"StretchBendAngles"
);
for
(
unsigned
int
ii
=
0
;
ii
<
(
int
)
bonds
.
getChildren
().
size
();
ii
++
)
{
for
(
auto
&
bond
:
bonds
.
getChildren
())
{
const
SerializationNode
&
bond
=
bonds
.
getChildren
()[
ii
];
double
k1
,
k2
;
double
k1
,
k2
;
if
(
version
==
1
)
if
(
version
==
1
)
k1
=
k2
=
bond
.
getDoubleProperty
(
"k"
);
k1
=
k2
=
bond
.
getDoubleProperty
(
"k"
);
...
...
plugins/cpupme/src/CpuPmeKernels.cpp
View file @
083bc501
...
@@ -515,8 +515,8 @@ CpuCalcPmeReciprocalForceKernel::~CpuCalcPmeReciprocalForceKernel() {
...
@@ -515,8 +515,8 @@ CpuCalcPmeReciprocalForceKernel::~CpuCalcPmeReciprocalForceKernel() {
pthread_mutex_destroy
(
&
lock
);
pthread_mutex_destroy
(
&
lock
);
pthread_cond_destroy
(
&
startCondition
);
pthread_cond_destroy
(
&
startCondition
);
pthread_cond_destroy
(
&
endCondition
);
pthread_cond_destroy
(
&
endCondition
);
for
(
int
i
=
0
;
i
<
(
int
)
tempGrid
.
size
();
i
++
)
for
(
auto
grid
:
tempGrid
)
fftwf_free
(
tempGrid
[
i
]
);
fftwf_free
(
grid
);
if
(
complexGrid
!=
NULL
)
if
(
complexGrid
!=
NULL
)
fftwf_free
(
complexGrid
);
fftwf_free
(
complexGrid
);
if
(
hasCreatedPlan
)
{
if
(
hasCreatedPlan
)
{
...
@@ -552,8 +552,8 @@ void CpuCalcPmeReciprocalForceKernel::runMainThread() {
...
@@ -552,8 +552,8 @@ void CpuCalcPmeReciprocalForceKernel::runMainThread() {
if
(
includeEnergy
)
{
if
(
includeEnergy
)
{
threads
.
resumeThreads
();
// Signal threads to compute energy.
threads
.
resumeThreads
();
// Signal threads to compute energy.
threads
.
waitForThreads
();
threads
.
waitForThreads
();
for
(
int
i
=
0
;
i
<
(
int
)
threadEnergy
.
size
();
i
++
)
for
(
auto
e
:
threadEnergy
)
energy
+=
threadEnergy
[
i
]
;
energy
+=
e
;
}
}
threads
.
resumeThreads
();
// Signal threads to perform reciprocal convolution.
threads
.
resumeThreads
();
// Signal threads to perform reciprocal convolution.
threads
.
waitForThreads
();
threads
.
waitForThreads
();
...
@@ -805,8 +805,8 @@ CpuCalcDispersionPmeReciprocalForceKernel::~CpuCalcDispersionPmeReciprocalForceK
...
@@ -805,8 +805,8 @@ CpuCalcDispersionPmeReciprocalForceKernel::~CpuCalcDispersionPmeReciprocalForceK
pthread_mutex_destroy
(
&
lock
);
pthread_mutex_destroy
(
&
lock
);
pthread_cond_destroy
(
&
startCondition
);
pthread_cond_destroy
(
&
startCondition
);
pthread_cond_destroy
(
&
endCondition
);
pthread_cond_destroy
(
&
endCondition
);
for
(
int
i
=
0
;
i
<
(
int
)
tempGrid
.
size
();
i
++
)
for
(
auto
grid
:
tempGrid
)
fftwf_free
(
tempGrid
[
i
]
);
fftwf_free
(
grid
);
if
(
complexGrid
!=
NULL
)
if
(
complexGrid
!=
NULL
)
fftwf_free
(
complexGrid
);
fftwf_free
(
complexGrid
);
if
(
hasCreatedPlan
)
{
if
(
hasCreatedPlan
)
{
...
@@ -843,8 +843,8 @@ void CpuCalcDispersionPmeReciprocalForceKernel::runMainThread() {
...
@@ -843,8 +843,8 @@ void CpuCalcDispersionPmeReciprocalForceKernel::runMainThread() {
if
(
includeEnergy
)
{
if
(
includeEnergy
)
{
threads
.
resumeThreads
();
// Signal threads to compute energy.
threads
.
resumeThreads
();
// Signal threads to compute energy.
threads
.
waitForThreads
();
threads
.
waitForThreads
();
for
(
int
i
=
0
;
i
<
(
int
)
threadEnergy
.
size
();
i
++
)
for
(
auto
e
:
threadEnergy
)
energy
+=
threadEnergy
[
i
]
;
energy
+=
e
;
}
}
threads
.
resumeThreads
();
// Signal threads to perform reciprocal convolution.
threads
.
resumeThreads
();
// Signal threads to perform reciprocal convolution.
threads
.
waitForThreads
();
threads
.
waitForThreads
();
...
...
plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp
View file @
083bc501
...
@@ -267,8 +267,7 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
...
@@ -267,8 +267,7 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
const
double
fscale
=
(
1
-
vscale
)
/
integrator
.
getFriction
();
const
double
fscale
=
(
1
-
vscale
)
/
integrator
.
getFriction
();
const
double
kT
=
BOLTZ
*
integrator
.
getTemperature
();
const
double
kT
=
BOLTZ
*
integrator
.
getTemperature
();
const
double
noisescale
=
sqrt
(
2
*
kT
*
integrator
.
getFriction
())
*
sqrt
(
0.5
*
(
1
-
vscale
*
vscale
)
/
integrator
.
getFriction
());
const
double
noisescale
=
sqrt
(
2
*
kT
*
integrator
.
getFriction
())
*
sqrt
(
0.5
*
(
1
-
vscale
*
vscale
)
/
integrator
.
getFriction
());
for
(
int
i
=
0
;
i
<
(
int
)
normalParticles
.
size
();
i
++
)
{
for
(
int
index
:
normalParticles
)
{
int
index
=
normalParticles
[
i
];
double
invMass
=
particleInvMass
[
index
];
double
invMass
=
particleInvMass
[
index
];
if
(
invMass
!=
0.0
)
{
if
(
invMass
!=
0.0
)
{
double
sqrtInvMass
=
sqrt
(
invMass
);
double
sqrtInvMass
=
sqrt
(
invMass
);
...
...
plugins/drude/serialization/src/DrudeForceProxy.cpp
View file @
083bc501
...
@@ -67,16 +67,12 @@ void* DrudeForceProxy::deserialize(const SerializationNode& node) const {
...
@@ -67,16 +67,12 @@ void* DrudeForceProxy::deserialize(const SerializationNode& node) const {
DrudeForce
*
force
=
new
DrudeForce
();
DrudeForce
*
force
=
new
DrudeForce
();
try
{
try
{
const
SerializationNode
&
particles
=
node
.
getChildNode
(
"Particles"
);
const
SerializationNode
&
particles
=
node
.
getChildNode
(
"Particles"
);
for
(
int
i
=
0
;
i
<
(
int
)
particles
.
getChildren
().
size
();
i
++
)
{
for
(
auto
&
particle
:
particles
.
getChildren
())
const
SerializationNode
&
particle
=
particles
.
getChildren
()[
i
];
force
->
addParticle
(
particle
.
getIntProperty
(
"p"
),
particle
.
getIntProperty
(
"p1"
),
particle
.
getIntProperty
(
"p2"
),
particle
.
getIntProperty
(
"p3"
),
particle
.
getIntProperty
(
"p4"
),
force
->
addParticle
(
particle
.
getIntProperty
(
"p"
),
particle
.
getIntProperty
(
"p1"
),
particle
.
getIntProperty
(
"p2"
),
particle
.
getIntProperty
(
"p3"
),
particle
.
getIntProperty
(
"p4"
),
particle
.
getDoubleProperty
(
"charge"
),
particle
.
getDoubleProperty
(
"polarizability"
),
particle
.
getDoubleProperty
(
"a12"
),
particle
.
getDoubleProperty
(
"a34"
));
particle
.
getDoubleProperty
(
"charge"
),
particle
.
getDoubleProperty
(
"polarizability"
),
particle
.
getDoubleProperty
(
"a12"
),
particle
.
getDoubleProperty
(
"a34"
));
}
const
SerializationNode
&
pairs
=
node
.
getChildNode
(
"ScreenedPairs"
);
const
SerializationNode
&
pairs
=
node
.
getChildNode
(
"ScreenedPairs"
);
for
(
int
i
=
0
;
i
<
(
int
)
pairs
.
getChildren
().
size
();
i
++
)
{
for
(
auto
&
pair
:
pairs
.
getChildren
())
const
SerializationNode
&
pair
=
pairs
.
getChildren
()[
i
];
force
->
addScreenedPair
(
pair
.
getIntProperty
(
"p1"
),
pair
.
getIntProperty
(
"p2"
),
pair
.
getDoubleProperty
(
"thole"
));
force
->
addScreenedPair
(
pair
.
getIntProperty
(
"p1"
),
pair
.
getIntProperty
(
"p2"
),
pair
.
getDoubleProperty
(
"thole"
));
}
}
}
catch
(...)
{
catch
(...)
{
delete
force
;
delete
force
;
...
...
plugins/rpmd/openmmapi/src/RPMDIntegrator.cpp
View file @
083bc501
...
@@ -119,13 +119,13 @@ State RPMDIntegrator::getState(int copy, int types, bool enforcePeriodicBox, int
...
@@ -119,13 +119,13 @@ State RPMDIntegrator::getState(int copy, int types, bool enforcePeriodicBox, int
const
vector
<
vector
<
int
>
>&
molecules
=
context
->
getMolecules
();
const
vector
<
vector
<
int
>
>&
molecules
=
context
->
getMolecules
();
Vec3
periodicBoxSize
[
3
];
Vec3
periodicBoxSize
[
3
];
state2
.
getPeriodicBoxVectors
(
periodicBoxSize
[
0
],
periodicBoxSize
[
1
],
periodicBoxSize
[
2
]);
state2
.
getPeriodicBoxVectors
(
periodicBoxSize
[
0
],
periodicBoxSize
[
1
],
periodicBoxSize
[
2
]);
for
(
int
i
=
0
;
i
<
(
int
)
molecules
.
size
();
i
++
)
{
for
(
auto
&
mol
:
molecules
)
{
// Find the molecule center.
// Find the molecule center.
Vec3
center
;
Vec3
center
;
for
(
int
j
=
0
;
j
<
(
int
)
molecules
[
i
].
size
();
j
++
)
for
(
int
j
:
mol
)
center
+=
refPos
[
molecules
[
i
][
j
]
];
center
+=
refPos
[
j
];
center
*=
1.0
/
mol
ecules
[
i
]
.
size
();
center
*=
1.0
/
mol
.
size
();
// Find the displacement to move it into the first periodic box.
// Find the displacement to move it into the first periodic box.
Vec3
diff
;
Vec3
diff
;
...
@@ -134,8 +134,8 @@ State RPMDIntegrator::getState(int copy, int types, bool enforcePeriodicBox, int
...
@@ -134,8 +134,8 @@ State RPMDIntegrator::getState(int copy, int types, bool enforcePeriodicBox, int
diff
+=
periodicBoxSize
[
0
]
*
floor
((
center
[
0
]
-
diff
[
0
])
/
periodicBoxSize
[
0
][
0
]);
diff
+=
periodicBoxSize
[
0
]
*
floor
((
center
[
0
]
-
diff
[
0
])
/
periodicBoxSize
[
0
][
0
]);
// Translate all the particles in the molecule.
// Translate all the particles in the molecule.
for
(
int
j
=
0
;
j
<
(
int
)
molecules
[
i
].
size
();
j
++
)
{
for
(
int
j
:
mol
)
{
Vec3
&
pos
=
positions
[
molecules
[
i
][
j
]
];
Vec3
&
pos
=
positions
[
j
];
pos
-=
diff
;
pos
-=
diff
;
}
}
}
}
...
@@ -188,9 +188,8 @@ void RPMDIntegrator::step(int steps) {
...
@@ -188,9 +188,8 @@ void RPMDIntegrator::step(int steps) {
context
->
getOwner
().
setPositions
(
p
);
context
->
getOwner
().
setPositions
(
p
);
isFirstStep
=
false
;
isFirstStep
=
false
;
}
}
vector
<
ForceImpl
*>&
forceImpls
=
context
->
getForceImpls
();
for
(
auto
impl
:
context
->
getForceImpls
())
{
for
(
int
i
=
0
;
i
<
(
int
)
forceImpls
.
size
();
i
++
)
{
RPMDUpdater
*
updater
=
dynamic_cast
<
RPMDUpdater
*>
(
impl
);
RPMDUpdater
*
updater
=
dynamic_cast
<
RPMDUpdater
*>
(
forceImpls
[
i
]);
if
(
updater
!=
NULL
)
if
(
updater
!=
NULL
)
updater
->
updateRPMDState
(
*
context
);
updater
->
updateRPMDState
(
*
context
);
}
}
...
...
plugins/rpmd/platforms/cuda/src/CudaRpmdKernels.cpp
View file @
083bc501
...
@@ -116,9 +116,9 @@ void CudaIntegrateRPMDStepKernel::initialize(const System& system, const RPMDInt
...
@@ -116,9 +116,9 @@ void CudaIntegrateRPMDStepKernel::initialize(const System& system, const RPMDInt
groupsNotContracted
=
-
1
;
groupsNotContracted
=
-
1
;
const
map
<
int
,
int
>&
contractions
=
integrator
.
getContractions
();
const
map
<
int
,
int
>&
contractions
=
integrator
.
getContractions
();
int
maxContractedCopies
=
0
;
int
maxContractedCopies
=
0
;
for
(
map
<
int
,
int
>::
const_iterator
iter
=
contractions
.
begin
();
iter
!=
contractions
.
end
();
++
iter
)
{
for
(
auto
&
c
:
contractions
)
{
int
group
=
iter
->
first
;
int
group
=
c
.
first
;
int
copies
=
iter
->
second
;
int
copies
=
c
.
second
;
if
(
group
<
0
||
group
>
31
)
if
(
group
<
0
||
group
>
31
)
throw
OpenMMException
(
"RPMDIntegrator: Force group must be between 0 and 31"
);
throw
OpenMMException
(
"RPMDIntegrator: Force group must be between 0 and 31"
);
if
(
copies
<
0
||
copies
>
numCopies
)
if
(
copies
<
0
||
copies
>
numCopies
)
...
@@ -166,8 +166,8 @@ void CudaIntegrateRPMDStepKernel::initialize(const System& system, const RPMDInt
...
@@ -166,8 +166,8 @@ void CudaIntegrateRPMDStepKernel::initialize(const System& system, const RPMDInt
// Create kernels for doing contractions.
// Create kernels for doing contractions.
for
(
map
<
int
,
int
>::
const_iterator
iter
=
groupsByCopies
.
begin
();
iter
!=
groupsByCopies
.
end
();
++
iter
)
{
for
(
auto
&
g
:
groupsByCopies
)
{
int
copies
=
iter
->
first
;
int
copies
=
g
.
first
;
replacements
.
clear
();
replacements
.
clear
();
replacements
[
"NUM_CONTRACTED_COPIES"
]
=
cu
.
intToString
(
copies
);
replacements
[
"NUM_CONTRACTED_COPIES"
]
=
cu
.
intToString
(
copies
);
replacements
[
"POS_SCALE"
]
=
cu
.
doubleToString
(
1.0
/
numCopies
);
replacements
[
"POS_SCALE"
]
=
cu
.
doubleToString
(
1.0
/
numCopies
);
...
@@ -267,9 +267,9 @@ void CudaIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
...
@@ -267,9 +267,9 @@ void CudaIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
// Now loop over contractions and compute forces from them.
// Now loop over contractions and compute forces from them.
for
(
map
<
int
,
int
>::
const_iterator
iter
=
groupsByCopies
.
begin
();
iter
!=
groupsByCopies
.
end
();
++
iter
)
{
for
(
auto
&
g
:
groupsByCopies
)
{
int
copies
=
iter
->
first
;
int
copies
=
g
.
first
;
int
groupFlags
=
iter
->
second
;
int
groupFlags
=
g
.
second
;
// Find the contracted positions.
// Find the contracted positions.
...
...
plugins/rpmd/platforms/opencl/src/OpenCLRpmdKernels.cpp
View file @
083bc501
...
@@ -96,9 +96,9 @@ void OpenCLIntegrateRPMDStepKernel::initialize(const System& system, const RPMDI
...
@@ -96,9 +96,9 @@ void OpenCLIntegrateRPMDStepKernel::initialize(const System& system, const RPMDI
groupsNotContracted
=
-
1
;
groupsNotContracted
=
-
1
;
const
map
<
int
,
int
>&
contractions
=
integrator
.
getContractions
();
const
map
<
int
,
int
>&
contractions
=
integrator
.
getContractions
();
int
maxContractedCopies
=
0
;
int
maxContractedCopies
=
0
;
for
(
map
<
int
,
int
>::
const_iterator
iter
=
contractions
.
begin
();
iter
!=
contractions
.
end
();
++
iter
)
{
for
(
auto
&
c
:
contractions
)
{
int
group
=
iter
->
first
;
int
group
=
c
.
first
;
int
copies
=
iter
->
second
;
int
copies
=
c
.
second
;
if
(
group
<
0
||
group
>
31
)
if
(
group
<
0
||
group
>
31
)
throw
OpenMMException
(
"RPMDIntegrator: Force group must be between 0 and 31"
);
throw
OpenMMException
(
"RPMDIntegrator: Force group must be between 0 and 31"
);
if
(
copies
<
0
||
copies
>
numCopies
)
if
(
copies
<
0
||
copies
>
numCopies
)
...
@@ -146,8 +146,8 @@ void OpenCLIntegrateRPMDStepKernel::initialize(const System& system, const RPMDI
...
@@ -146,8 +146,8 @@ void OpenCLIntegrateRPMDStepKernel::initialize(const System& system, const RPMDI
// Create kernels for doing contractions.
// Create kernels for doing contractions.
for
(
map
<
int
,
int
>::
const_iterator
iter
=
groupsByCopies
.
begin
();
iter
!=
groupsByCopies
.
end
();
++
iter
)
{
for
(
auto
&
g
:
groupsByCopies
)
{
int
copies
=
iter
->
first
;
int
copies
=
g
.
first
;
replacements
.
clear
();
replacements
.
clear
();
replacements
[
"NUM_CONTRACTED_COPIES"
]
=
cl
.
intToString
(
copies
);
replacements
[
"NUM_CONTRACTED_COPIES"
]
=
cl
.
intToString
(
copies
);
replacements
[
"POS_SCALE"
]
=
cl
.
doubleToString
(
1.0
/
numCopies
);
replacements
[
"POS_SCALE"
]
=
cl
.
doubleToString
(
1.0
/
numCopies
);
...
@@ -182,8 +182,8 @@ void OpenCLIntegrateRPMDStepKernel::initializeKernels(ContextImpl& context) {
...
@@ -182,8 +182,8 @@ void OpenCLIntegrateRPMDStepKernel::initializeKernels(ContextImpl& context) {
copyFromContextKernel
.
setArg
<
cl
::
Buffer
>
(
3
,
velocities
->
getDeviceBuffer
());
copyFromContextKernel
.
setArg
<
cl
::
Buffer
>
(
3
,
velocities
->
getDeviceBuffer
());
copyFromContextKernel
.
setArg
<
cl
::
Buffer
>
(
4
,
cl
.
getPosq
().
getDeviceBuffer
());
copyFromContextKernel
.
setArg
<
cl
::
Buffer
>
(
4
,
cl
.
getPosq
().
getDeviceBuffer
());
copyFromContextKernel
.
setArg
<
cl
::
Buffer
>
(
6
,
cl
.
getAtomIndexArray
().
getDeviceBuffer
());
copyFromContextKernel
.
setArg
<
cl
::
Buffer
>
(
6
,
cl
.
getAtomIndexArray
().
getDeviceBuffer
());
for
(
map
<
int
,
int
>::
const_iterator
iter
=
groupsByCopies
.
begin
();
iter
!=
groupsByCopies
.
end
();
++
iter
)
{
for
(
auto
&
g
:
groupsByCopies
)
{
int
copies
=
iter
->
first
;
int
copies
=
g
.
first
;
positionContractionKernels
[
copies
].
setArg
<
cl
::
Buffer
>
(
0
,
positions
->
getDeviceBuffer
());
positionContractionKernels
[
copies
].
setArg
<
cl
::
Buffer
>
(
0
,
positions
->
getDeviceBuffer
());
positionContractionKernels
[
copies
].
setArg
<
cl
::
Buffer
>
(
1
,
contractedPositions
->
getDeviceBuffer
());
positionContractionKernels
[
copies
].
setArg
<
cl
::
Buffer
>
(
1
,
contractedPositions
->
getDeviceBuffer
());
forceContractionKernels
[
copies
].
setArg
<
cl
::
Buffer
>
(
0
,
forces
->
getDeviceBuffer
());
forceContractionKernels
[
copies
].
setArg
<
cl
::
Buffer
>
(
0
,
forces
->
getDeviceBuffer
());
...
@@ -286,9 +286,9 @@ void OpenCLIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
...
@@ -286,9 +286,9 @@ void OpenCLIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
copyToContextKernel
.
setArg
<
cl
::
Buffer
>
(
2
,
contractedPositions
->
getDeviceBuffer
());
copyToContextKernel
.
setArg
<
cl
::
Buffer
>
(
2
,
contractedPositions
->
getDeviceBuffer
());
copyFromContextKernel
.
setArg
<
cl
::
Buffer
>
(
1
,
contractedForces
->
getDeviceBuffer
());
copyFromContextKernel
.
setArg
<
cl
::
Buffer
>
(
1
,
contractedForces
->
getDeviceBuffer
());
copyFromContextKernel
.
setArg
<
cl
::
Buffer
>
(
5
,
contractedPositions
->
getDeviceBuffer
());
copyFromContextKernel
.
setArg
<
cl
::
Buffer
>
(
5
,
contractedPositions
->
getDeviceBuffer
());
for
(
map
<
int
,
int
>::
const_iterator
iter
=
groupsByCopies
.
begin
();
iter
!=
groupsByCopies
.
end
();
++
iter
)
{
for
(
auto
&
g
:
groupsByCopies
)
{
int
copies
=
iter
->
first
;
int
copies
=
g
.
first
;
int
groupFlags
=
iter
->
second
;
int
groupFlags
=
g
.
second
;
// Find the contracted positions.
// Find the contracted positions.
...
...
plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.cpp
View file @
083bc501
...
@@ -55,9 +55,9 @@ static vector<Vec3>& extractForces(ContextImpl& context) {
...
@@ -55,9 +55,9 @@ static vector<Vec3>& extractForces(ContextImpl& context) {
ReferenceIntegrateRPMDStepKernel
::~
ReferenceIntegrateRPMDStepKernel
()
{
ReferenceIntegrateRPMDStepKernel
::~
ReferenceIntegrateRPMDStepKernel
()
{
if
(
fft
!=
NULL
)
if
(
fft
!=
NULL
)
fftpack_destroy
(
fft
);
fftpack_destroy
(
fft
);
for
(
map
<
int
,
fftpack
*>::
const_iterator
iter
=
contractionFFT
.
begin
();
iter
!=
contractionFFT
.
end
();
++
iter
)
for
(
auto
&
c
:
contractionFFT
)
if
(
iter
->
second
!=
NULL
)
if
(
c
.
second
!=
NULL
)
fftpack_destroy
(
iter
->
second
);
fftpack_destroy
(
c
.
second
);
}
}
void
ReferenceIntegrateRPMDStepKernel
::
initialize
(
const
System
&
system
,
const
RPMDIntegrator
&
integrator
)
{
void
ReferenceIntegrateRPMDStepKernel
::
initialize
(
const
System
&
system
,
const
RPMDIntegrator
&
integrator
)
{
...
@@ -79,9 +79,9 @@ void ReferenceIntegrateRPMDStepKernel::initialize(const System& system, const RP
...
@@ -79,9 +79,9 @@ void ReferenceIntegrateRPMDStepKernel::initialize(const System& system, const RP
groupsNotContracted
=
-
1
;
groupsNotContracted
=
-
1
;
const
map
<
int
,
int
>&
contractions
=
integrator
.
getContractions
();
const
map
<
int
,
int
>&
contractions
=
integrator
.
getContractions
();
int
maxContractedCopies
=
0
;
int
maxContractedCopies
=
0
;
for
(
map
<
int
,
int
>::
const_iterator
iter
=
contractions
.
begin
();
iter
!=
contractions
.
end
();
++
iter
)
{
for
(
auto
&
c
:
contractions
)
{
int
group
=
iter
->
first
;
int
group
=
c
.
first
;
int
copies
=
iter
->
second
;
int
copies
=
c
.
second
;
if
(
group
<
0
||
group
>
31
)
if
(
group
<
0
||
group
>
31
)
throw
OpenMMException
(
"RPMDIntegrator: Force group must be between 0 and 31"
);
throw
OpenMMException
(
"RPMDIntegrator: Force group must be between 0 and 31"
);
if
(
copies
<
0
||
copies
>
numCopies
)
if
(
copies
<
0
||
copies
>
numCopies
)
...
@@ -290,9 +290,9 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const
...
@@ -290,9 +290,9 @@ void ReferenceIntegrateRPMDStepKernel::computeForces(ContextImpl& context, const
// Now loop over contractions and compute forces from them.
// Now loop over contractions and compute forces from them.
for
(
map
<
int
,
int
>::
const_iterator
iter
=
groupsByCopies
.
begin
();
iter
!=
groupsByCopies
.
end
();
++
iter
)
{
for
(
auto
&
g
:
groupsByCopies
)
{
int
copies
=
iter
->
first
;
int
copies
=
g
.
first
;
int
groupFlags
=
iter
->
second
;
int
groupFlags
=
g
.
second
;
fftpack
*
shortFFT
=
contractionFFT
[
copies
];
fftpack
*
shortFFT
=
contractionFFT
[
copies
];
// Find the contracted positions.
// Find the contracted positions.
...
...
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