NonbondedForceImpl.cpp 16.4 KB
Newer Older
1
2
3
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
Evan Pretti's avatar
Evan Pretti committed
4
5
 * This is part of the OpenMM molecular simulation toolkit.                   *
 * See https://openmm.org/development.                                        *
6
 *                                                                            *
7
 * Portions copyright (c) 2008-2024 Stanford University and the Authors.      *
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 * Authors: Peter Eastman                                                     *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

30
31
32
#ifdef WIN32
  #define _USE_MATH_DEFINES // Needed to get M_PI
#endif
33
#include "openmm/OpenMMException.h"
34
#include "openmm/internal/ContextImpl.h"
35
#include "openmm/internal/Messages.h"
36
37
#include "openmm/internal/NonbondedForceImpl.h"
#include "openmm/kernels.h"
38
#include <cmath>
39
#include <map>
40
#include <sstream>
41
#include <algorithm>
42
43

using namespace OpenMM;
44
using namespace std;
45

46
NonbondedForceImpl::NonbondedForceImpl(const NonbondedForce& owner) : owner(owner) {
47
48
49
50
51
    forceGroup = owner.getForceGroup();
    recipForceGroup = owner.getReciprocalSpaceForceGroup();
    if (recipForceGroup < 0)
        recipForceGroup = owner.getForceGroup();
    includeDirectSpace = owner.getIncludeDirectSpace();
52
53
}

54
NonbondedForceImpl::~NonbondedForceImpl() {
55
56
}

57
void NonbondedForceImpl::initialize(ContextImpl& context) {
58
    kernel = context.getPlatform().createKernel(CalcNonbondedForceKernel::Name(), context);
Peter Eastman's avatar
Peter Eastman committed
59

60
    // Check for errors in the specification of exceptions.
Peter Eastman's avatar
Peter Eastman committed
61

62
    const System& system = context.getSystem();
63
64
    if (owner.getNumParticles() != system.getNumParticles())
        throw OpenMMException("NonbondedForce must have exactly as many particles as the System it belongs to.");
65
66
67
68
    if (owner.getUseSwitchingFunction()) {
        if (owner.getSwitchingDistance() < 0 || owner.getSwitchingDistance() >= owner.getCutoffDistance())
            throw OpenMMException("NonbondedForce: Switching distance must satisfy 0 <= r_switch < r_cutoff");
    }
69
70
71
72
73
74
75
76
    for (int i = 0; i < owner.getNumParticles(); i++) {
        double charge, sigma, epsilon;
        owner.getParticleParameters(i, charge, sigma, epsilon);
        if (sigma < 0)
            throw OpenMMException("NonbondedForce: sigma for a particle cannot be negative");
        if (epsilon < 0)
            throw OpenMMException("NonbondedForce: epsilon for a particle cannot be negative");
    }
77
78
    vector<set<int> > exceptions(owner.getNumParticles());
    for (int i = 0; i < owner.getNumExceptions(); i++) {
79
        int particle[2];
80
        double chargeProd, sigma, epsilon;
81
        owner.getExceptionParameters(i, particle[0], particle[1], chargeProd, sigma, epsilon);
82
83
        int minp = min(particle[0], particle[1]);
        int maxp = max(particle[0], particle[1]);
84
85
86
87
88
89
90
91
        for (int j = 0; j < 2; j++) {
            if (particle[j] < 0 || particle[j] >= owner.getNumParticles()) {
                stringstream msg;
                msg << "NonbondedForce: Illegal particle index for an exception: ";
                msg << particle[j];
                throw OpenMMException(msg.str());
            }
        }
92
        if (exceptions[minp].count(maxp) > 0) {
93
            stringstream msg;
94
95
96
97
            msg << "NonbondedForce: Multiple exceptions are specified for particles ";
            msg << particle[0];
            msg << " and ";
            msg << particle[1];
98
99
            throw OpenMMException(msg.str());
        }
100
        exceptions[minp].insert(maxp);
101
102
        if (sigma < 0)
            throw OpenMMException("NonbondedForce: sigma for an exception cannot be negative");
103
104
105
106
107
108
109
110
111
        if (epsilon < 0)
            throw OpenMMException("NonbondedForce: epsilon for an exception cannot be negative");
    }
    for (int i = 0; i < owner.getNumParticleParameterOffsets(); i++) {
        string parameter;
        int particleIndex;
        double chargeScale, sigmaScale, epsilonScale;
        owner.getParticleParameterOffset(i, parameter, particleIndex, chargeScale, sigmaScale, epsilonScale);
        if (particleIndex < 0 || particleIndex >= owner.getNumParticles()) {
112
            stringstream msg;
113
114
            msg << "NonbondedForce: Illegal particle index for a particle parameter offset: ";
            msg << particleIndex;
115
116
            throw OpenMMException(msg.str());
        }
117
118
119
120
121
122
123
    }
    for (int i = 0; i < owner.getNumExceptionParameterOffsets(); i++) {
        string parameter;
        int exceptionIndex;
        double chargeScale, sigmaScale, epsilonScale;
        owner.getExceptionParameterOffset(i, parameter, exceptionIndex, chargeScale, sigmaScale, epsilonScale);
        if (exceptionIndex < 0 || exceptionIndex >= owner.getNumExceptions()) {
124
            stringstream msg;
125
126
            msg << "NonbondedForce: Illegal exception index for an exception parameter offset: ";
            msg << exceptionIndex;
127
128
            throw OpenMMException(msg.str());
        }
Peter Eastman's avatar
Peter Eastman committed
129
    }
Peter Eastman's avatar
Peter Eastman committed
130
    if (owner.getNonbondedMethod() != NonbondedForce::NoCutoff && owner.getNonbondedMethod() != NonbondedForce::CutoffNonPeriodic) {
131
        Vec3 boxVectors[3];
132
        system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
133
134
        double cutoff = owner.getCutoffDistance();
        if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
135
            throw OpenMMException("NonbondedForce: "+Messages::cutoffTooLarge);
136
137
        if (owner.getNonbondedMethod() == NonbondedForce::Ewald && (boxVectors[1][0] != 0.0 || boxVectors[2][0] != 0.0 || boxVectors[2][1] != 0))
            throw OpenMMException("NonbondedForce: Ewald is not supported with non-rectangular boxes.  Use PME instead.");
138
    }
139
    kernel.getAs<CalcNonbondedForceKernel>().initialize(context.getSystem(), owner);
140
141
}

142
double NonbondedForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
143
144
    bool includeDirect = (includeDirectSpace && (groups&(1<<forceGroup)) != 0);
    bool includeReciprocal = ((groups&(1<<recipForceGroup)) != 0);
145
    return kernel.getAs<CalcNonbondedForceKernel>().execute(context, includeForces, includeEnergy, includeDirect, includeReciprocal);
146
147
}

148
map<string, double> NonbondedForceImpl::getDefaultParameters() {
149
150
151
152
    map<string, double> parameters;
    for (int i = 0; i < owner.getNumGlobalParameters(); i++)
        parameters[owner.getGlobalParameterName(i)] = owner.getGlobalParameterDefaultValue(i);
    return parameters;
153
154
}

155
std::vector<std::string> NonbondedForceImpl::getKernelNames() {
156
    std::vector<std::string> names;
157
    names.push_back(CalcNonbondedForceKernel::Name());
158
159
    return names;
}
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179

class NonbondedForceImpl::ErrorFunction {
public:
    virtual double getValue(int arg) const = 0;
};

class NonbondedForceImpl::EwaldErrorFunction : public ErrorFunction {
public:
    EwaldErrorFunction(double width, double alpha, double target) : width(width), alpha(alpha), target(target) {
    }
    double getValue(int arg) const {
        double temp = arg*M_PI/(width*alpha);
        return target-0.05*sqrt(width*alpha)*arg*exp(-temp*temp);
    }
private:
    double width, alpha, target;
};

void NonbondedForceImpl::calcEwaldParameters(const System& system, const NonbondedForce& force, double& alpha, int& kmaxx, int& kmaxy, int& kmaxz) {
    Vec3 boxVectors[3];
180
    system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
181
    double tol = force.getEwaldErrorTolerance();
Peter Eastman's avatar
Peter Eastman committed
182
    alpha = (1.0/force.getCutoffDistance())*std::sqrt(-log(2.0*tol));
183
184
185
186
187
188
189
190
191
192
193
    kmaxx = findZero(EwaldErrorFunction(boxVectors[0][0], alpha, tol), 10);
    kmaxy = findZero(EwaldErrorFunction(boxVectors[1][1], alpha, tol), 10);
    kmaxz = findZero(EwaldErrorFunction(boxVectors[2][2], alpha, tol), 10);
    if (kmaxx%2 == 0)
        kmaxx++;
    if (kmaxy%2 == 0)
        kmaxy++;
    if (kmaxz%2 == 0)
        kmaxz++;
}

194
195
void NonbondedForceImpl::calcPMEParameters(const System& system, const NonbondedForce& force, double& alpha, int& xsize, int& ysize, int& zsize, bool lj) {
    if (lj)
196
        force.getLJPMEParameters(alpha, xsize, ysize, zsize);
197
    else
198
        force.getPMEParameters(alpha, xsize, ysize, zsize);
199
200
201
202
203
    if (alpha == 0.0) {
        Vec3 boxVectors[3];
        system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
        double tol = force.getEwaldErrorTolerance();
        alpha = (1.0/force.getCutoffDistance())*std::sqrt(-log(2.0*tol));
Peter Eastman's avatar
Peter Eastman committed
204
        if (lj) {
205
206
207
            xsize = (int) ceil(alpha*boxVectors[0][0]/(3*pow(tol, 0.2)));
            ysize = (int) ceil(alpha*boxVectors[1][1]/(3*pow(tol, 0.2)));
            zsize = (int) ceil(alpha*boxVectors[2][2]/(3*pow(tol, 0.2)));
Peter Eastman's avatar
Peter Eastman committed
208
209
210
211
212
213
214
215
216
        }
        else {
            xsize = (int) ceil(2*alpha*boxVectors[0][0]/(3*pow(tol, 0.2)));
            ysize = (int) ceil(2*alpha*boxVectors[1][1]/(3*pow(tol, 0.2)));
            zsize = (int) ceil(2*alpha*boxVectors[2][2]/(3*pow(tol, 0.2)));
        }
        xsize = max(xsize, 6);
        ysize = max(ysize, 6);
        zsize = max(zsize, 6);
217
    }
218
219
}

220
int NonbondedForceImpl::findZero(const NonbondedForceImpl::ErrorFunction& f, int initialGuess) {
221
222
223
224
225
226
227
228
229
230
231
    int arg = initialGuess;
    double value = f.getValue(arg);
    if (value > 0.0) {
        while (value > 0.0 && arg > 0)
            value = f.getValue(--arg);
        return arg+1;
    }
    while (value < 0.0)
        value = f.getValue(++arg);
    return arg;
}
232

233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
double NonbondedForceImpl::evalIntegral(double r, double rs, double rc, double sigma) {
    // Compute the indefinite integral of the LJ interaction multiplied by the switching function.
    // This is a large and somewhat horrifying expression, though it does grow on you if you look
    // at it long enough.  Perhaps it could be simplified further, but I got tired of working on it.
    
    double A = 1/(rc-rs);
    double A2 = A*A;
    double A3 = A2*A;
    double sig2 = sigma*sigma;
    double sig6 = sig2*sig2*sig2;
    double rs2 = rs*rs;
    double rs3 = rs*rs2;
    double r2 = r*r;
    double r3 = r*r2;
    double r4 = r*r3;
    double r5 = r*r4;
    double r6 = r*r5;
    double r9 = r3*r6;
    return sig6*A3*((
        sig6*(
            + rs3*28*(6*rs2*A2 + 15*rs*A + 10)
            - r*rs2*945*(rs2*A2 + 2*rs*A + 1)
            + r2*rs*1080*(2*rs2*A2 + 3*rs*A + 1)
            - r3*420*(6*rs2*A2 + 6*rs*A + 1)
            + r4*756*(2*rs*A2 + A)
            - r5*378*A2)
        -r6*(
            + rs3*84*(6*rs2*A2 + 15*rs*A + 10)
            - r*rs2*3780*(rs2*A2 + 2*rs*A + 1)
            + r2*rs*7560*(2*rs2*A2 + 3*rs*A + 1))
        )/(252*r9)
     - log(r)*10*(6*rs2*A2 + 6*rs*A + 1)
     + r*15*(2*rs*A2 + A)
     - r2*3*A2
    );
}

270
double NonbondedForceImpl::calcDispersionCorrection(const System& system, const NonbondedForce& force) {
271
272
    if (force.getNonbondedMethod() == NonbondedForce::NoCutoff || force.getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
        return 0.0;
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293

    // 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;
    }

294
295
296
297
298
    // Identify all particle classes (defined by sigma and epsilon), and count the number of
    // particles in each class.

    map<pair<double, double>, int> classCounts;
    for (int i = 0; i < force.getNumParticles(); i++) {
299
        pair<double, double> key = make_pair(sigma[i], epsilon[i]);
300
301
302
303
304
305
306
307
308
        map<pair<double, double>, int>::iterator entry = classCounts.find(key);
        if (entry == classCounts.end())
            classCounts[key] = 1;
        else
            entry->second++;
    }

    // Loop over all pairs of classes to compute the coefficient.

309
310
311
312
    double sum1 = 0, sum2 = 0, sum3 = 0;
    bool useSwitch = force.getUseSwitchingFunction();
    double cutoff = force.getCutoffDistance();
    double switchDist = force.getSwitchingDistance();
313
314
315
    for (map<pair<double, double>, int>::const_iterator entry = classCounts.begin(); entry != classCounts.end(); ++entry) {
        double sigma = entry->first.first;
        double epsilon = entry->first.second;
316
317
        double count = (double) entry->second;
        count *= (count + 1) / 2;
318
319
320
321
        double sigma2 = sigma*sigma;
        double sigma6 = sigma2*sigma2*sigma2;
        sum1 += count*epsilon*sigma6*sigma6;
        sum2 += count*epsilon*sigma6;
322
323
        if (useSwitch)
            sum3 += count*epsilon*(evalIntegral(cutoff, switchDist, cutoff, sigma)-evalIntegral(switchDist, switchDist, cutoff, sigma));
324
325
326
327
328
    }
    for (map<pair<double, double>, int>::const_iterator class1 = classCounts.begin(); class1 != classCounts.end(); ++class1)
        for (map<pair<double, double>, int>::const_iterator class2 = classCounts.begin(); class2 != class1; ++class2) {
            double sigma = 0.5*(class1->first.first+class2->first.first);
            double epsilon = sqrt(class1->first.second*class2->first.second);
329
330
            double count = (double) class1->second;
            count *= (double) class2->second;
331
332
333
334
            double sigma2 = sigma*sigma;
            double sigma6 = sigma2*sigma2*sigma2;
            sum1 += count*epsilon*sigma6*sigma6;
            sum2 += count*epsilon*sigma6;
335
336
            if (useSwitch)
                sum3 += count*epsilon*(evalIntegral(cutoff, switchDist, cutoff, sigma)-evalIntegral(switchDist, switchDist, cutoff, sigma));
337
        }
338
    double numParticles = (double) system.getNumParticles();
Jason Swails's avatar
Jason Swails committed
339
    double numInteractions = (numParticles*(numParticles+1))/2;
340
341
    sum1 /= numInteractions;
    sum2 /= numInteractions;
342
343
    sum3 /= numInteractions;
    return 8*numParticles*numParticles*M_PI*(sum1/(9*pow(cutoff, 9))-sum2/(3*pow(cutoff, 3))+sum3);
344
}
345

346
347
void NonbondedForceImpl::updateParametersInContext(ContextImpl& context, int firstParticle, int lastParticle, int firstException, int lastException) {
    kernel.getAs<CalcNonbondedForceKernel>().copyParametersToContext(context, owner, firstParticle, lastParticle, firstException, lastException);
348
    context.systemChanged();
349
}
350
351
352
353

void NonbondedForceImpl::getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
    kernel.getAs<CalcNonbondedForceKernel>().getPMEParameters(alpha, nx, ny, nz);
}
354
355
356
357

void NonbondedForceImpl::getLJPMEParameters(double& alpha, int& nx, int& ny, int& nz) const {
    kernel.getAs<CalcNonbondedForceKernel>().getLJPMEParameters(alpha, nx, ny, nz);
}