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
030ef272
Commit
030ef272
authored
Jan 13, 2010
by
Mark Friedrichs
Browse files
Remove unused methods
parent
074f6b36
Changes
2
Show whitespace changes
Inline
Side-by-side
Showing
2 changed files
with
0 additions
and
697 deletions
+0
-697
platforms/reference/src/gbsa/CpuObc.cpp
platforms/reference/src/gbsa/CpuObc.cpp
+0
-639
platforms/reference/src/gbsa/CpuObc.h
platforms/reference/src/gbsa/CpuObc.h
+0
-58
No files found.
platforms/reference/src/gbsa/CpuObc.cpp
View file @
030ef272
...
@@ -674,642 +674,3 @@ std::string CpuObc::getStateString( const char* title ) const {
...
@@ -674,642 +674,3 @@ std::string CpuObc::getStateString( const char* title ) const {
return
message
.
str
();
return
message
.
str
();
}
}
/**---------------------------------------------------------------------------------------
Write Born energy and forces (Simbios)
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
@param resultsFileName output file name
@return SimTKOpenMMCommon::DefaultReturn unless
file cannot be opened
in which case return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int
CpuObc
::
writeBornEnergyForces
(
RealOpenMM
**
atomCoordinates
,
const
RealOpenMM
*
partialCharges
,
RealOpenMM
**
forces
,
const
std
::
string
&
resultsFileName
)
const
{
// ---------------------------------------------------------------------------------------
static
const
char
*
methodName
=
"
\n
CpuObc::writeBornEnergyForces"
;
// ---------------------------------------------------------------------------------------
ImplicitSolventParameters
*
implicitSolventParameters
=
getImplicitSolventParameters
();
const
ObcParameters
*
obcParameters
=
static_cast
<
const
ObcParameters
*>
(
implicitSolventParameters
);
int
numberOfAtoms
=
obcParameters
->
getNumberOfAtoms
();
const
RealOpenMM
*
atomicRadii
=
obcParameters
->
getAtomicRadii
();
const
RealOpenMM
*
bornRadii
=
getBornRadiiConst
();
const
RealOpenMM
*
scaledRadii
=
obcParameters
->
getScaledRadiusFactors
();
const
RealOpenMM
*
obcChain
=
getObcChainConst
();
const
RealOpenMM
energy
=
getEnergy
();
// ---------------------------------------------------------------------------------------
// open file -- return if unsuccessful
FILE
*
implicitSolventResultsFile
=
NULL
;
#ifdef _MSC_VER
fopen_s
(
&
implicitSolventResultsFile
,
resultsFileName
.
c_str
(),
"w"
);
#else
implicitSolventResultsFile
=
fopen
(
resultsFileName
.
c_str
(),
"w"
);
#endif
// diganostics
std
::
stringstream
message
;
message
<<
methodName
;
if
(
implicitSolventResultsFile
!=
NULL
){
std
::
stringstream
message
;
message
<<
methodName
;
message
<<
" Opened file=<"
<<
resultsFileName
<<
">."
;
SimTKOpenMMLog
::
printMessage
(
message
);
}
else
{
std
::
stringstream
message
;
message
<<
methodName
;
message
<<
" could not open file=<"
<<
resultsFileName
<<
"> -- abort output."
;
SimTKOpenMMLog
::
printMessage
(
message
);
return
SimTKOpenMMCommon
::
ErrorReturn
;
}
// header
(
void
)
fprintf
(
implicitSolventResultsFile
,
"# %d atoms E=%.7e format: coords(3) bornRadii(input) q atomicRadii scaleFactors forces obcChain
\n
"
,
numberOfAtoms
,
energy
);
RealOpenMM
forceConversion
=
(
RealOpenMM
)
1.0
;
RealOpenMM
lengthConversion
=
(
RealOpenMM
)
1.0
;
// output
if
(
forces
!=
NULL
&&
atomCoordinates
!=
NULL
&&
partialCharges
!=
NULL
&&
atomicRadii
!=
NULL
){
for
(
int
ii
=
0
;
ii
<
numberOfAtoms
;
ii
++
){
(
void
)
fprintf
(
implicitSolventResultsFile
,
"%.7e %.7e %.7e %.7e %.5f %.5f %.5f %.7e %.7e %.7e %.7e
\n
"
,
lengthConversion
*
atomCoordinates
[
ii
][
0
],
lengthConversion
*
atomCoordinates
[
ii
][
1
],
lengthConversion
*
atomCoordinates
[
ii
][
2
],
(
bornRadii
!=
NULL
?
lengthConversion
*
bornRadii
[
ii
]
:
0.0
),
partialCharges
[
ii
],
lengthConversion
*
atomicRadii
[
ii
],
scaledRadii
[
ii
],
forceConversion
*
forces
[
ii
][
0
],
forceConversion
*
forces
[
ii
][
1
],
forceConversion
*
forces
[
ii
][
2
],
forceConversion
*
obcChain
[
ii
]
);
}
}
(
void
)
fclose
(
implicitSolventResultsFile
);
return
SimTKOpenMMCommon
::
DefaultReturn
;
}
/**---------------------------------------------------------------------------------------
Write results from first loop
@param numberOfAtoms number of atoms
@param forces forces
@param bornForce Born force prefactor
@param outputFileName output file name
@return SimTKOpenMMCommon::DefaultReturn unless
file cannot be opened
in which case return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int
CpuObc
::
writeForceLoop1
(
int
numberOfAtoms
,
RealOpenMM
**
forces
,
const
RealOpenMM
*
bornForce
,
const
std
::
string
&
outputFileName
){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::writeForceLoop1";
// ---------------------------------------------------------------------------------------
int
chunkSize
;
if
(
bornForce
){
chunkSize
=
3
;
}
else
{
chunkSize
=
4
;
}
StringVector
lineVector
;
std
::
stringstream
header
;
lineVector
.
push_back
(
"# bornF F"
);
for
(
int
atomI
=
0
;
atomI
<
numberOfAtoms
;
atomI
++
){
std
::
stringstream
line
;
line
<<
(
atomI
+
1
)
<<
" "
;
SimTKOpenMMUtilities
::
formatRealStringStream
(
line
,
forces
[
atomI
],
chunkSize
);
if
(
bornForce
){
line
<<
" "
<<
bornForce
[
atomI
];
}
lineVector
.
push_back
(
line
.
str
()
);
}
return
SimTKOpenMMUtilities
::
writeFile
(
lineVector
,
outputFileName
);
}
/**---------------------------------------------------------------------------------------
Write results
@param numberOfAtoms number of atoms
@param chunkSizes vector of chunk sizes for realRealOpenMMVector
@param realRealOpenMMVector vector of RealOpenMM**
@param realVector vector of RealOpenMM*
@param outputFileName output file name
@return SimTKOpenMMCommon::DefaultReturn unless
file cannot be opened
in which case return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int
CpuObc
::
writeForceLoop
(
int
numberOfAtoms
,
const
IntVector
&
chunkSizes
,
const
RealOpenMMPtrPtrVector
&
realRealOpenMMVector
,
const
RealOpenMMPtrVector
&
realVector
,
const
std
::
string
&
outputFileName
){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::writeForceLoop";
static
const
int
maxChunks
=
10
;
int
chunks
[
maxChunks
];
// ---------------------------------------------------------------------------------------
for
(
int
ii
=
0
;
ii
<
(
int
)
chunkSizes
.
size
();
ii
++
){
chunks
[
ii
]
=
chunkSizes
[
ii
];
}
for
(
int
ii
=
(
int
)
chunkSizes
.
size
();
ii
<
maxChunks
;
ii
++
){
chunks
[
ii
]
=
3
;
}
StringVector
lineVector
;
std
::
stringstream
header
;
// lineVector.push_back( "# " );
for
(
int
atomI
=
0
;
atomI
<
numberOfAtoms
;
atomI
++
){
std
::
stringstream
line
;
char
buffer
[
128
];
(
void
)
sprintf
(
buffer
,
"%4d "
,
atomI
);
line
<<
buffer
;
int
index
=
0
;
for
(
RealOpenMMPtrPtrVectorCI
ii
=
realRealOpenMMVector
.
begin
();
ii
!=
realRealOpenMMVector
.
end
();
ii
++
){
RealOpenMM
**
forces
=
*
ii
;
(
void
)
sprintf
(
buffer
,
"%11.5f %11.5f %11.5f "
,
forces
[
atomI
][
0
],
forces
[
atomI
][
1
],
forces
[
atomI
][
2
]
);
line
<<
buffer
;
// SimTKOpenMMUtilities::formatRealStringStream( line, forces[atomI], chunks[index++] );
// line << " ";
}
for
(
RealOpenMMPtrVectorCI
ii
=
realVector
.
begin
();
ii
!=
realVector
.
end
();
ii
++
){
RealOpenMM
*
array
=
*
ii
;
(
void
)
sprintf
(
buffer
,
"%11.5f "
,
array
[
atomI
]
);
line
<<
buffer
;
}
lineVector
.
push_back
(
line
.
str
()
);
}
return
SimTKOpenMMUtilities
::
writeFile
(
lineVector
,
outputFileName
);
}
/**---------------------------------------------------------------------------------------
Get Obc Born energy and forces -- used debugging
@param bornRadii Born radii -- optional; if NULL, then ObcParameters
entry is used
@param atomCoordinates atomic coordinates
@param partialCharges partial charges
@param forces forces
@return SimTKOpenMMCommon::DefaultReturn;
The array bornRadii is also updated and the obcEnergy
--------------------------------------------------------------------------------------- */
int
CpuObc
::
computeBornEnergyForcesPrint
(
RealOpenMM
*
bornRadii
,
RealOpenMM
**
atomCoordinates
,
const
RealOpenMM
*
partialCharges
,
RealOpenMM
**
forces
){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nCpuObc::computeBornEnergyForcesPrint";
static
const
RealOpenMM
zero
=
(
RealOpenMM
)
0.0
;
static
const
RealOpenMM
one
=
(
RealOpenMM
)
1.0
;
static
const
RealOpenMM
two
=
(
RealOpenMM
)
2.0
;
static
const
RealOpenMM
three
=
(
RealOpenMM
)
3.0
;
static
const
RealOpenMM
four
=
(
RealOpenMM
)
4.0
;
static
const
RealOpenMM
half
=
(
RealOpenMM
)
0.5
;
static
const
RealOpenMM
fourth
=
(
RealOpenMM
)
0.25
;
static
const
RealOpenMM
eighth
=
(
RealOpenMM
)
0.125
;
// ---------------------------------------------------------------------------------------
const
ObcParameters
*
obcParameters
=
getObcParameters
();
const
int
numberOfAtoms
=
obcParameters
->
getNumberOfAtoms
();
if
(
bornRadii
==
NULL
){
bornRadii
=
getBornRadii
();
}
// suppress warning about fopen in Visual Studio
#if defined(_MSC_VER)
#pragma warning(push)
#pragma warning(disable:4996)
#endif
FILE
*
logFile
=
NULL
;
//FILE* logFile = SimTKOpenMMLog::getSimTKOpenMMLogFile( );
//FILE* logFile = fopen( "bF", "w" );
#if defined(_MSC_VER)
#pragma warning(pop)
#endif
// ---------------------------------------------------------------------------------------
// constants
const
RealOpenMM
preFactor
=
obcParameters
->
getPreFactor
();
const
RealOpenMM
dielectricOffset
=
obcParameters
->
getDielectricOffset
();
// ---------------------------------------------------------------------------------------
// set energy/forces to zero
RealOpenMM
obcEnergy
=
zero
;
const
unsigned
int
arraySzInBytes
=
sizeof
(
RealOpenMM
)
*
numberOfAtoms
;
for
(
int
ii
=
0
;
ii
<
numberOfAtoms
;
ii
++
){
memset
(
forces
[
ii
],
0
,
3
*
sizeof
(
RealOpenMM
)
);
}
RealOpenMM
*
bornForces
=
getBornForce
();
memset
(
bornForces
,
0
,
arraySzInBytes
);
// ---------------------------------------------------------------------------------------
// N*( 8 + pow) ACE
// compute the nonpolar solvation via ACE approximation
if
(
includeAceApproximation
()
){
computeAceNonPolarForce
(
obcParameters
,
bornRadii
,
&
obcEnergy
,
bornForces
);
if
(
logFile
){
(
void
)
fprintf
(
logFile
,
"
\n
ACE E=%.5e
\n
"
,
obcEnergy
);
for
(
int
atomI
=
0
;
atomI
<
numberOfAtoms
;
atomI
++
){
(
void
)
fprintf
(
logFile
,
" %d bR=%.6e bF=%.6e
\n
"
,
atomI
,
bornRadii
[
atomI
],
bornForces
[
atomI
]
);
}
}
}
// ---------------------------------------------------------------------------------------
// first main loop
for
(
int
atomI
=
0
;
atomI
<
numberOfAtoms
;
atomI
++
){
RealOpenMM
partialChargeI
=
preFactor
*
partialCharges
[
atomI
];
for
(
int
atomJ
=
atomI
;
atomJ
<
numberOfAtoms
;
atomJ
++
){
RealOpenMM
deltaR
[
ReferenceForce
::
LastDeltaRIndex
];
if
(
_obcParameters
->
getPeriodic
())
ReferenceForce
::
getDeltaRPeriodic
(
atomCoordinates
[
atomI
],
atomCoordinates
[
atomJ
],
_obcParameters
->
getPeriodicBox
(),
deltaR
);
else
ReferenceForce
::
getDeltaR
(
atomCoordinates
[
atomI
],
atomCoordinates
[
atomJ
],
deltaR
);
if
(
_obcParameters
->
getUseCutoff
()
&&
deltaR
[
ReferenceForce
::
RIndex
]
>
_obcParameters
->
getCutoffDistance
())
continue
;
RealOpenMM
r2
=
deltaR
[
ReferenceForce
::
R2Index
];
RealOpenMM
deltaX
=
deltaR
[
ReferenceForce
::
XIndex
];
RealOpenMM
deltaY
=
deltaR
[
ReferenceForce
::
YIndex
];
RealOpenMM
deltaZ
=
deltaR
[
ReferenceForce
::
ZIndex
];
// 3 FLOP
RealOpenMM
alpha2_ij
=
bornRadii
[
atomI
]
*
bornRadii
[
atomJ
];
RealOpenMM
D_ij
=
r2
/
(
four
*
alpha2_ij
);
// exp + 2 + sqrt FLOP
RealOpenMM
expTerm
=
EXP
(
-
D_ij
);
RealOpenMM
denominator2
=
r2
+
alpha2_ij
*
expTerm
;
RealOpenMM
denominator
=
SQRT
(
denominator2
);
// 6 FLOP
RealOpenMM
Gpol
=
(
partialChargeI
*
partialCharges
[
atomJ
])
/
denominator
;
// dGpol/dr = -1/2*(Gpol/denominator2)*(2r - r/2*exp() )
RealOpenMM
dGpol_dr
=
-
Gpol
*
(
one
-
fourth
*
expTerm
)
/
denominator2
;
// 5 FLOP
RealOpenMM
dGpol_dalpha2_ij
=
-
half
*
Gpol
*
expTerm
*
(
one
+
D_ij
)
/
denominator2
;
// 11 FLOP
if
(
atomI
!=
atomJ
){
bornForces
[
atomJ
]
+=
dGpol_dalpha2_ij
*
bornRadii
[
atomI
];
deltaX
*=
dGpol_dr
;
deltaY
*=
dGpol_dr
;
deltaZ
*=
dGpol_dr
;
forces
[
atomI
][
0
]
+=
deltaX
;
forces
[
atomI
][
1
]
+=
deltaY
;
forces
[
atomI
][
2
]
+=
deltaZ
;
forces
[
atomJ
][
0
]
-=
deltaX
;
forces
[
atomJ
][
1
]
-=
deltaY
;
forces
[
atomJ
][
2
]
-=
deltaZ
;
}
else
{
Gpol
*=
half
;
}
// 3 FLOP
obcEnergy
+=
Gpol
;
bornForces
[
atomI
]
+=
dGpol_dalpha2_ij
*
bornRadii
[
atomJ
];
//if( logFile && (atomI == -1 || atomJ == -1) ){
// (void) fprintf( logFile, "\nWWX %d %d F[%.6e %.6e %.6e] bF=[%.6e %.6e] Gpl[%.6e %.6e %.6e] rb[%6.4f %7.4f] rs[%6.4f %7.4f] ",
// atomI, atomJ,
// forces[atomI][0], forces[atomI][1], forces[atomI][2],
// bornForces[atomI], bornForces[atomJ],
// Gpol,dGpol_dr,dGpol_dalpha2_ij,
// bornRadii[atomI],bornRadii[atomJ],atomicRadii[atomI],atomicRadii[atomJ] );
//
// (void) fprintf( logFile, "\nWWX %d %d %.1f r2=%.4f q=%.2f bF=[%.6e %.6e] Gpl[%.6e %.6e %.6e] rb[%.5f %.5f] add[%.6e %.6e] ",
// atomI, atomJ, preFactor, r2, partialCharges[atomJ],
// bornForces[atomI], bornForces[atomJ],
// Gpol,dGpol_dr,dGpol_dalpha2_ij,
// bornRadii[atomI], bornRadii[atomJ],
// dGpol_dalpha2_ij*bornRadii[atomJ], dGpol_dalpha2_ij*bornRadii[atomI] );
//}
}
}
if
(
logFile
){
(
void
)
fprintf
(
logFile
,
"
\n
WXX bF & F E=%.8e preFactor=%.5f"
,
obcEnergy
,
preFactor
);
for
(
int
atomI
=
0
;
atomI
<
numberOfAtoms
;
atomI
++
){
(
void
)
fprintf
(
logFile
,
"
\n
WXX %d q=%.4f bR=%.5e bF=%.3f F[%.6e %.6e %.6e] "
,
atomI
,
partialCharges
[
atomI
],
bornRadii
[
atomI
],
bornForces
[
atomI
],
forces
[
atomI
][
0
],
forces
[
atomI
][
1
],
forces
[
atomI
][
2
]
);
}
}
if
(
1
){
std
::
string
outputFileName
=
"Loop1Cpu.txt"
;
CpuObc
::
writeForceLoop1
(
numberOfAtoms
,
forces
,
bornForces
,
outputFileName
);
/*
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
forces[atomI][0] = forces[atomI][1] = forces[atomI][2] = (RealOpenMM) 0.0;
}
*/
}
// ---------------------------------------------------------------------------------------
// second main loop
// initialize Born radii & ObcChain temp arrays -- contain values
// used in next iteration
RealOpenMM
*
bornRadiiTemp
=
getBornRadiiTemp
();
memset
(
bornRadiiTemp
,
0
,
arraySzInBytes
);
RealOpenMM
*
obcChainTemp
=
getObcChainTemp
();
memset
(
obcChainTemp
,
0
,
arraySzInBytes
);
RealOpenMM
*
obcChain
=
getObcChain
();
const
RealOpenMM
*
atomicRadii
=
obcParameters
->
getAtomicRadii
();
const
RealOpenMM
alphaObc
=
obcParameters
->
getAlphaObc
();
const
RealOpenMM
betaObc
=
obcParameters
->
getBetaObc
();
const
RealOpenMM
gammaObc
=
obcParameters
->
getGammaObc
();
const
RealOpenMM
*
scaledRadiusFactor
=
obcParameters
->
getScaledRadiusFactors
();
// compute factor that depends only on the outer loop index
for
(
int
atomI
=
0
;
atomI
<
numberOfAtoms
;
atomI
++
){
bornForces
[
atomI
]
*=
bornRadii
[
atomI
]
*
bornRadii
[
atomI
]
*
obcChain
[
atomI
];
}
if
(
1
){
std
::
string
outputFileName
=
"PostLoop1Cpu.txt"
;
IntVector
chunkVector
;
chunkVector
.
push_back
(
3
);
RealOpenMMPtrPtrVector
realPtrPtrVector
;
realPtrPtrVector
.
push_back
(
forces
);
RealOpenMMPtrVector
realPtrVector
;
realPtrVector
.
push_back
(
bornRadii
);
realPtrVector
.
push_back
(
bornForces
);
realPtrVector
.
push_back
(
obcChain
);
CpuObc
::
writeForceLoop
(
numberOfAtoms
,
chunkVector
,
realPtrPtrVector
,
realPtrVector
,
outputFileName
);
}
RealOpenMM
*
bornSumArray
=
(
RealOpenMM
*
)
malloc
(
sizeof
(
RealOpenMM
)
*
numberOfAtoms
);
memset
(
bornSumArray
,
0
,
sizeof
(
RealOpenMM
)
*
numberOfAtoms
);
/*
for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
forces[atomI][0] = 0.0;
forces[atomI][1] = 0.0;
forces[atomI][2] = 0.0;
} */
for
(
int
atomI
=
0
;
atomI
<
numberOfAtoms
;
atomI
++
){
// radius w/ dielectric offset applied
RealOpenMM
radiusI
=
atomicRadii
[
atomI
];
RealOpenMM
offsetRadiusI
=
radiusI
-
dielectricOffset
;
// used to compute Born radius for next iteration
RealOpenMM
bornSum
=
zero
;
for
(
int
atomJ
=
0
;
atomJ
<
numberOfAtoms
;
atomJ
++
){
if
(
atomJ
!=
atomI
){
RealOpenMM
deltaR
[
ReferenceForce
::
LastDeltaRIndex
];
if
(
_obcParameters
->
getPeriodic
())
ReferenceForce
::
getDeltaRPeriodic
(
atomCoordinates
[
atomI
],
atomCoordinates
[
atomJ
],
_obcParameters
->
getPeriodicBox
(),
deltaR
);
else
ReferenceForce
::
getDeltaR
(
atomCoordinates
[
atomI
],
atomCoordinates
[
atomJ
],
deltaR
);
if
(
_obcParameters
->
getUseCutoff
()
&&
deltaR
[
ReferenceForce
::
RIndex
]
>
_obcParameters
->
getCutoffDistance
())
continue
;
RealOpenMM
deltaX
=
deltaR
[
ReferenceForce
::
XIndex
];
RealOpenMM
deltaY
=
deltaR
[
ReferenceForce
::
YIndex
];
RealOpenMM
deltaZ
=
deltaR
[
ReferenceForce
::
ZIndex
];
RealOpenMM
r
=
deltaR
[
ReferenceForce
::
RIndex
];
// radius w/ dielectric offset applied
RealOpenMM
radiusJ
=
atomicRadii
[
atomJ
]
-
dielectricOffset
;
RealOpenMM
scaledRadiusJ
=
radiusJ
*
scaledRadiusFactor
[
atomJ
];
RealOpenMM
scaledRadiusJ2
=
scaledRadiusJ
*
scaledRadiusJ
;
RealOpenMM
rScaledRadiusJ
=
r
+
scaledRadiusJ
;
// L_ij != 1 && U_ij != 1
// dL/dr & dU/dr are zero (this can be shown analytically)
// removed from calculation
if
(
offsetRadiusI
<
rScaledRadiusJ
){
RealOpenMM
l_ij
=
offsetRadiusI
>
FABS
(
r
-
scaledRadiusJ
)
?
offsetRadiusI
:
FABS
(
r
-
scaledRadiusJ
);
l_ij
=
one
/
l_ij
;
RealOpenMM
l_ij2
=
l_ij
*
l_ij
;
RealOpenMM
u_ij
=
one
/
rScaledRadiusJ
;
RealOpenMM
u_ij2
=
u_ij
*
u_ij
;
RealOpenMM
rInverse
=
one
/
r
;
RealOpenMM
r2Inverse
=
rInverse
*
rInverse
;
RealOpenMM
logRatio
=
LN
(
u_ij
/
l_ij
);
RealOpenMM
t3
=
eighth
*
(
one
+
scaledRadiusJ2
*
r2Inverse
)
*
(
l_ij2
-
u_ij2
)
+
fourth
*
logRatio
*
r2Inverse
;
RealOpenMM
de
=
bornForces
[
atomI
]
*
t3
*
rInverse
;
deltaX
*=
de
;
deltaY
*=
de
;
deltaZ
*=
de
;
forces
[
atomI
][
0
]
-=
deltaX
;
forces
[
atomI
][
1
]
-=
deltaY
;
forces
[
atomI
][
2
]
-=
deltaZ
;
forces
[
atomJ
][
0
]
+=
deltaX
;
forces
[
atomJ
][
1
]
+=
deltaY
;
forces
[
atomJ
][
2
]
+=
deltaZ
;
// Born radius term
RealOpenMM
term
=
l_ij
-
u_ij
+
fourth
*
r
*
(
u_ij2
-
l_ij2
)
+
(
half
*
rInverse
)
*
logRatio
+
(
fourth
*
scaledRadiusJ
*
scaledRadiusJ
*
rInverse
)
*
(
l_ij2
-
u_ij2
);
if
(
offsetRadiusI
<
(
scaledRadiusJ
-
r
)
){
term
+=
two
*
(
(
one
/
offsetRadiusI
)
-
l_ij
);
}
bornSum
+=
term
;
if
(
atomI
==
-
1
||
atomJ
==
-
1
){
(
void
)
fprintf
(
logFile
,
"
\n
XXY %d %d de=%.6e bF[%.6e %6e] t3=%.6e r=%.6e trm=%.6e bSm=%.6e f[%.6e %.6e %.6e]"
,
atomI
,
atomJ
,
de
,
bornForces
[
atomI
],
obcChain
[
atomI
],
t3
,
r
,
term
,
bornSum
,
forces
[
atomI
][
0
],
forces
[
atomI
][
1
],
forces
[
atomI
][
2
]
);
}
}
}
}
bornSumArray
[
atomI
]
=
bornSum
;
// OBC-specific code (Eqs. 6-8 in paper)
bornSum
*=
half
*
offsetRadiusI
;
RealOpenMM
sum2
=
bornSum
*
bornSum
;
RealOpenMM
sum3
=
bornSum
*
sum2
;
RealOpenMM
tanhSum
=
TANH
(
alphaObc
*
bornSum
-
betaObc
*
sum2
+
gammaObc
*
sum3
);
bornRadiiTemp
[
atomI
]
=
one
/
(
one
/
offsetRadiusI
-
tanhSum
/
radiusI
);
obcChainTemp
[
atomI
]
=
offsetRadiusI
*
(
alphaObc
-
two
*
betaObc
*
bornSum
+
three
*
gammaObc
*
sum2
);
obcChainTemp
[
atomI
]
=
(
one
-
tanhSum
*
tanhSum
)
*
obcChainTemp
[
atomI
]
/
radiusI
;
if
(
logFile
&&
atomI
>=
0
){
(
void
)
fprintf
(
logFile
,
"
\n
XXX %d bSum[%.6e %.6e %.6e] bRt=[%.6e %6e] obc=%.6e rI=[%.5f %.5f]"
,
atomI
,
bornSumArray
[
atomI
],
bornSum
,
tanhSum
,
bornRadii
[
atomI
],
bornRadiiTemp
[
atomI
],
obcChainTemp
[
atomI
],
radiusI
,
offsetRadiusI
);
}
}
RealOpenMM
conversion
=
(
RealOpenMM
)
1
;
for
(
int
atomI
=
0
;
atomI
<
numberOfAtoms
;
atomI
++
){
forces
[
atomI
][
0
]
*=
conversion
;
forces
[
atomI
][
1
]
*=
conversion
;
forces
[
atomI
][
2
]
*=
conversion
;
}
setEnergy
(
obcEnergy
*
conversion
);
if
(
1
){
std
::
string
outputFileName
=
"Loop2Cpu.txt"
;
IntVector
chunkVector
;
chunkVector
.
push_back
(
3
);
RealOpenMMPtrPtrVector
realPtrPtrVector
;
realPtrPtrVector
.
push_back
(
forces
);
RealOpenMMPtrVector
realPtrVector
;
realPtrVector
.
push_back
(
bornSumArray
);
// realPtrVector.push_back( bornRadiiTemp );
// realPtrVector.push_back( obcChainTemp );
CpuObc
::
writeForceLoop
(
numberOfAtoms
,
chunkVector
,
realPtrPtrVector
,
realPtrVector
,
outputFileName
);
}
if
(
bornSumArray
){
free
(
(
char
*
)
bornSumArray
);
}
// 6 FLOP
/*
RealOpenMM forceFactor = getForceConversionFactor();
RealOpenMM constantFactor = 1.0f/electricConstant;
if( fabs(forceFactor - 1.0f) > 1.0e-04 ){
constantFactor *= forceFactor;
for( int ii = 0; ii < numberOfAtoms; ii++ ){
forces[ii][0] *= forceFactor;
forces[ii][1] *= forceFactor;
forces[ii][2] *= forceFactor;
}
} */
// copy new Born radii and obcChain values into permanent array
//(void) fprintf( logFile, "\nBorn radii not being updated!!!!" );
memcpy
(
bornRadii
,
bornRadiiTemp
,
arraySzInBytes
);
memcpy
(
obcChain
,
obcChainTemp
,
arraySzInBytes
);
if
(
logFile
){
(
void
)
fclose
(
logFile
);
}
return
SimTKOpenMMCommon
::
DefaultReturn
;
}
platforms/reference/src/gbsa/CpuObc.h
View file @
030ef272
...
@@ -147,9 +147,6 @@ class CpuObc : public CpuImplicitSolvent {
...
@@ -147,9 +147,6 @@ class CpuObc : public CpuImplicitSolvent {
int
computeBornEnergyForces
(
RealOpenMM
*
bornRadii
,
RealOpenMM
**
atomCoordinates
,
int
computeBornEnergyForces
(
RealOpenMM
*
bornRadii
,
RealOpenMM
**
atomCoordinates
,
const
RealOpenMM
*
partialCharges
,
RealOpenMM
**
forces
);
const
RealOpenMM
*
partialCharges
,
RealOpenMM
**
forces
);
int
computeBornEnergyForcesPrint
(
RealOpenMM
*
bornRadii
,
RealOpenMM
**
atomCoordinates
,
const
RealOpenMM
*
partialCharges
,
RealOpenMM
**
forces
);
/**---------------------------------------------------------------------------------------
/**---------------------------------------------------------------------------------------
Get state
Get state
...
@@ -162,61 +159,6 @@ class CpuObc : public CpuImplicitSolvent {
...
@@ -162,61 +159,6 @@ class CpuObc : public CpuImplicitSolvent {
std
::
string
getStateString
(
const
char
*
title
)
const
;
std
::
string
getStateString
(
const
char
*
title
)
const
;
/**---------------------------------------------------------------------------------------
Write Born energy and forces (Simbios)
@param atomCoordinates atomic coordinates
@param partialCharges partial atom charges
@param forces force array
@param resultsFileName output file name
@return SimTKOpenMMCommon::DefaultReturn if file opened; else return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
int
writeBornEnergyForces
(
RealOpenMM
**
atomCoordinates
,
const
RealOpenMM
*
partialCharges
,
RealOpenMM
**
forces
,
const
std
::
string
&
resultsFileName
)
const
;
/**---------------------------------------------------------------------------------------
Write results from first loop
@param atomCoordinates atomic coordinates
@param RealOpenMM forces forces
@param outputFileName output file name
@return SimTKOpenMMCommon::DefaultReturn unless
file cannot be opened
in which case return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
static
int
writeForceLoop1
(
int
numberOfAtoms
,
RealOpenMM
**
forces
,
const
RealOpenMM
*
bornForce
,
const
std
::
string
&
outputFileName
);
/**---------------------------------------------------------------------------------------
Write results
@param numberOfAtoms number of atoms
@param chunkSizes vector of chunk sizes for realRealOpenMMVector
@param realRealOpenMMVector vector of RealOpenMM**
@param realVector vector of RealOpenMM*
@param outputFileName output file name
@return SimTKOpenMMCommon::DefaultReturn unless
file cannot be opened
in which case return SimTKOpenMMCommon::ErrorReturn
--------------------------------------------------------------------------------------- */
static
int
writeForceLoop
(
int
numberOfAtoms
,
const
IntVector
&
chunkSizes
,
const
RealOpenMMPtrPtrVector
&
realRealOpenMMVector
,
const
RealOpenMMPtrVector
&
realVector
,
const
std
::
string
&
outputFileName
);
};
};
// ---------------------------------------------------------------------------------------
// ---------------------------------------------------------------------------------------
...
...
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