Commit 61b3cd18 authored by peastman's avatar peastman
Browse files

Dispersion correction is based on default values of offset parameters

parent aed9eed8
...@@ -236,14 +236,32 @@ double NonbondedForceImpl::calcDispersionCorrection(const System& system, const ...@@ -236,14 +236,32 @@ double NonbondedForceImpl::calcDispersionCorrection(const System& system, const
if (force.getNonbondedMethod() == NonbondedForce::NoCutoff || force.getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic) if (force.getNonbondedMethod() == NonbondedForce::NoCutoff || force.getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
return 0.0; return 0.0;
// Record sigma and epsilon for every particle, including the default value
// for every offset parameter.
vector<double> sigma(force.getNumParticles()), epsilon(force.getNumParticles());
for (int i = 0; i < force.getNumParticles(); i++) {
double charge;
force.getParticleParameters(i, charge, sigma[i], epsilon[i]);
}
map<string, double> param;
for (int i = 0; i < force.getNumGlobalParameters(); i++)
param[force.getGlobalParameterName(i)] = force.getGlobalParameterDefaultValue(i);
for (int i = 0; i < force.getNumParticleParameterOffsets(); i++) {
string parameter;
int index;
double chargeScale, sigmaScale, epsilonScale;
force.getParticleParameterOffset(i, parameter, index, chargeScale, sigmaScale, epsilonScale);
sigma[index] += param[parameter]*sigmaScale;
epsilon[index] += param[parameter]*epsilonScale;
}
// Identify all particle classes (defined by sigma and epsilon), and count the number of // Identify all particle classes (defined by sigma and epsilon), and count the number of
// particles in each class. // particles in each class.
map<pair<double, double>, int> classCounts; map<pair<double, double>, int> classCounts;
for (int i = 0; i < force.getNumParticles(); i++) { for (int i = 0; i < force.getNumParticles(); i++) {
double charge, sigma, epsilon; pair<double, double> key = make_pair(sigma[i], epsilon[i]);
force.getParticleParameters(i, charge, sigma, epsilon);
pair<double, double> key = make_pair(sigma, epsilon);
map<pair<double, double>, int>::iterator entry = classCounts.find(key); map<pair<double, double>, int>::iterator entry = classCounts.find(key);
if (entry == classCounts.end()) if (entry == classCounts.end())
classCounts[key] = 1; classCounts[key] = 1;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment