CustomNonbondedForce.cpp 15 KB
Newer Older
1
2
3
4
5
6
7
8
/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
9
 * Portions copyright (c) 2008-2024 Stanford University and the Authors.      *
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
 * 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.                                     *
 * -------------------------------------------------------------------------- */

#include "openmm/Force.h"
#include "openmm/OpenMMException.h"
#include "openmm/CustomNonbondedForce.h"
35
#include "openmm/internal/AssertionUtilities.h"
36
37
38
39
40
41
42
#include "openmm/internal/CustomNonbondedForceImpl.h"
#include <cmath>
#include <map>
#include <sstream>
#include <utility>

using namespace OpenMM;
43
using namespace std;
44

45
CustomNonbondedForce::CustomNonbondedForce(const string& energy) : energyExpression(energy), nonbondedMethod(NoCutoff), cutoffDistance(1.0),
46
    switchingDistance(-1.0), useSwitchingFunction(false), useLongRangeCorrection(false), numContexts(0) {
47
48
}

49
CustomNonbondedForce::CustomNonbondedForce(const CustomNonbondedForce& rhs) : numContexts(0) {
50
    // Copy everything and deep copy the tabulated functions
51
52
    setForceGroup(rhs.getForceGroup());
    setName(rhs.getName());
53
54
55
56
57
58
59
60
    energyExpression = rhs.energyExpression;
    nonbondedMethod = rhs.nonbondedMethod;
    cutoffDistance = rhs.cutoffDistance;
    switchingDistance = rhs.switchingDistance;
    useSwitchingFunction = rhs.useSwitchingFunction;
    useLongRangeCorrection = rhs.useLongRangeCorrection;
    parameters = rhs.parameters;
    globalParameters = rhs.globalParameters;
61
    computedValues = rhs.computedValues;
62
    energyParameterDerivatives = rhs.energyParameterDerivatives;
63
64
65
    particles = rhs.particles;
    exclusions = rhs.exclusions;
    interactionGroups = rhs.interactionGroups;
66
    for (vector<FunctionInfo>::const_iterator it = rhs.functions.begin(); it != rhs.functions.end(); ++it)
67
        functions.push_back(FunctionInfo(it->name, it->function->Copy()));
68
69
}

70
CustomNonbondedForce::~CustomNonbondedForce() {
peastman's avatar
peastman committed
71
72
    for (auto function : functions)
        delete function.function;
73
74
}

75
76
77
78
79
80
81
82
83
84
85
86
87
const string& CustomNonbondedForce::getEnergyFunction() const {
    return energyExpression;
}

void CustomNonbondedForce::setEnergyFunction(const std::string& energy) {
    energyExpression = energy;
}

CustomNonbondedForce::NonbondedMethod CustomNonbondedForce::getNonbondedMethod() const {
    return nonbondedMethod;
}

void CustomNonbondedForce::setNonbondedMethod(NonbondedMethod method) {
88
89
    if (method < 0 || method > 2)
        throw OpenMMException("CustomNonbondedForce: Illegal value for nonbonded method");
90
91
92
93
94
95
96
97
98
99
100
    nonbondedMethod = method;
}

double CustomNonbondedForce::getCutoffDistance() const {
    return cutoffDistance;
}

void CustomNonbondedForce::setCutoffDistance(double distance) {
    cutoffDistance = distance;
}

101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
bool CustomNonbondedForce::getUseSwitchingFunction() const {
    return useSwitchingFunction;
}

void CustomNonbondedForce::setUseSwitchingFunction(bool use) {
    useSwitchingFunction = use;
}

double CustomNonbondedForce::getSwitchingDistance() const {
    return switchingDistance;
}

void CustomNonbondedForce::setSwitchingDistance(double distance) {
    switchingDistance = distance;
}

bool CustomNonbondedForce::getUseLongRangeCorrection() const {
    return useLongRangeCorrection;
}

void CustomNonbondedForce::setUseLongRangeCorrection(bool use) {
    useLongRangeCorrection = use;
}

125
126
int CustomNonbondedForce::addPerParticleParameter(const string& name) {
    parameters.push_back(PerParticleParameterInfo(name));
127
128
129
    return parameters.size()-1;
}

130
const string& CustomNonbondedForce::getPerParticleParameterName(int index) const {
131
    ASSERT_VALID_INDEX(index, parameters);
132
133
134
    return parameters[index].name;
}

135
void CustomNonbondedForce::setPerParticleParameterName(int index, const string& name) {
136
    ASSERT_VALID_INDEX(index, parameters);
137
    parameters[index].name = name;
138
139
}

140
141
int CustomNonbondedForce::addGlobalParameter(const string& name, double defaultValue) {
    globalParameters.push_back(GlobalParameterInfo(name, defaultValue));
142
143
144
145
    return globalParameters.size()-1;
}

const string& CustomNonbondedForce::getGlobalParameterName(int index) const {
146
    ASSERT_VALID_INDEX(index, globalParameters);
147
    return globalParameters[index].name;
148
149
150
}

void CustomNonbondedForce::setGlobalParameterName(int index, const string& name) {
151
    ASSERT_VALID_INDEX(index, globalParameters);
152
153
154
155
    globalParameters[index].name = name;
}

double CustomNonbondedForce::getGlobalParameterDefaultValue(int index) const {
156
    ASSERT_VALID_INDEX(index, globalParameters);
157
158
159
160
    return globalParameters[index].defaultValue;
}

void CustomNonbondedForce::setGlobalParameterDefaultValue(int index, double defaultValue) {
161
    ASSERT_VALID_INDEX(index, globalParameters);
162
    globalParameters[index].defaultValue = defaultValue;
163
164
}

165
166
167
168
169
170
171
172
173
174
175
176
177
178
void CustomNonbondedForce::addEnergyParameterDerivative(const string& name) {
    for (int i = 0; i < globalParameters.size(); i++)
        if (name == globalParameters[i].name) {
            energyParameterDerivatives.push_back(i);
            return;
        }
    throw OpenMMException(string("addEnergyParameterDerivative: Unknown global parameter '"+name+"'"));
}

const string& CustomNonbondedForce::getEnergyParameterDerivativeName(int index) const {
    ASSERT_VALID_INDEX(index, energyParameterDerivatives);
    return globalParameters[energyParameterDerivatives[index]].name;
}

179
180
181
182
183
184
int CustomNonbondedForce::addParticle(const vector<double>& parameters) {
    particles.push_back(ParticleInfo(parameters));
    return particles.size()-1;
}

void CustomNonbondedForce::getParticleParameters(int index, std::vector<double>& parameters) const {
185
    ASSERT_VALID_INDEX(index, particles);
186
187
188
189
    parameters = particles[index].parameters;
}

void CustomNonbondedForce::setParticleParameters(int index, const vector<double>& parameters) {
190
    ASSERT_VALID_INDEX(index, particles);
191
    particles[index].parameters = parameters;
192
193
194
195
    if (numContexts > 0) {
        firstChangedParticle = min(index, firstChangedParticle);
        lastChangedParticle = max(index, lastChangedParticle);
    }
196
197
}

198
199
200
201
202
int CustomNonbondedForce::addExclusion(int particle1, int particle2) {
    exclusions.push_back(ExclusionInfo(particle1, particle2));
    return exclusions.size()-1;
}
void CustomNonbondedForce::getExclusionParticles(int index, int& particle1, int& particle2) const {
203
    ASSERT_VALID_INDEX(index, exclusions);
204
205
206
207
208
    particle1 = exclusions[index].particle1;
    particle2 = exclusions[index].particle2;
}

void CustomNonbondedForce::setExclusionParticles(int index, int particle1, int particle2) {
209
    ASSERT_VALID_INDEX(index, exclusions);
210
211
    exclusions[index].particle1 = particle1;
    exclusions[index].particle2 = particle2;
212
}
213
214

void CustomNonbondedForce::createExclusionsFromBonds(const vector<pair<int, int> >& bonds, int bondCutoff) {
215
216
    if (bondCutoff < 1)
        return;
peastman's avatar
peastman committed
217
218
    for (auto& bond : bonds)
        if (bond.first < 0 || bond.second < 0 || bond.first >= particles.size() || bond.second >= particles.size())
219
            throw OpenMMException("createExclusionsFromBonds: Illegal particle index in list of bonds");
220
221
    vector<set<int> > exclusions(particles.size());
    vector<set<int> > bonded12(exclusions.size());
peastman's avatar
peastman committed
222
223
224
    for (auto& bond : bonds) {
        int p1 = bond.first;
        int p2 = bond.second;
225
226
227
228
229
230
231
        exclusions[p1].insert(p2);
        exclusions[p2].insert(p1);
        bonded12[p1].insert(p2);
        bonded12[p2].insert(p1);
    }
    for (int level = 0; level < bondCutoff-1; level++) {
        vector<set<int> > currentExclusions = exclusions;
peastman's avatar
peastman committed
232
233
234
        for (int i = 0; i < (int) particles.size(); i++)
            for (int j : currentExclusions[i])
                exclusions[j].insert(bonded12[i].begin(), bonded12[i].end());
235
236
    }
    for (int i = 0; i < (int) exclusions.size(); ++i)
peastman's avatar
peastman committed
237
238
239
        for (int j : exclusions[i])
            if (j < i)
                addExclusion(j, i);
240
241
}

242
int CustomNonbondedForce::addTabulatedFunction(const std::string& name, TabulatedFunction* function) {
243
244
245
246
    functions.push_back(FunctionInfo(name, function));
    return functions.size()-1;
}

247
const TabulatedFunction& CustomNonbondedForce::getTabulatedFunction(int index) const {
248
249
250
251
    ASSERT_VALID_INDEX(index, functions);
    return *functions[index].function;
}

252
TabulatedFunction& CustomNonbondedForce::getTabulatedFunction(int index) {
253
254
255
256
    ASSERT_VALID_INDEX(index, functions);
    return *functions[index].function;
}

257
const string& CustomNonbondedForce::getTabulatedFunctionName(int index) const {
258
259
260
    ASSERT_VALID_INDEX(index, functions);
    return functions[index].name;
}
261

262
int CustomNonbondedForce::addFunction(const std::string& name, const std::vector<double>& values, double min, double max) {
263
    functions.push_back(FunctionInfo(name, new Continuous1DFunction(values, min, max)));
264
265
266
    return functions.size()-1;
}

267
void CustomNonbondedForce::getFunctionParameters(int index, std::string& name, std::vector<double>& values, double& min, double& max) const {
268
    ASSERT_VALID_INDEX(index, functions);
269
270
271
    Continuous1DFunction* function = dynamic_cast<Continuous1DFunction*>(functions[index].function);
    if (function == NULL)
        throw OpenMMException("CustomNonbondedForce: function is not a Continuous1DFunction");
272
    name = functions[index].name;
273
    function->getFunctionParameters(values, min, max);
274
275
}

276
void CustomNonbondedForce::setFunctionParameters(int index, const std::string& name, const std::vector<double>& values, double min, double max) {
277
    ASSERT_VALID_INDEX(index, functions);
278
279
280
    Continuous1DFunction* function = dynamic_cast<Continuous1DFunction*>(functions[index].function);
    if (function == NULL)
        throw OpenMMException("CustomNonbondedForce: function is not a Continuous1DFunction");
281
    functions[index].name = name;
282
    function->setFunctionParameters(values, min, max);
283
284
}

285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
int CustomNonbondedForce::addComputedValue(const string& name, const string& expression) {
    computedValues.push_back(ComputedValueInfo(name, expression));
    return computedValues.size()-1;
}

void CustomNonbondedForce::getComputedValueParameters(int index, string& name, string& expression) const {
    ASSERT_VALID_INDEX(index, computedValues);
    name = computedValues[index].name;
    expression = computedValues[index].expression;
}

void CustomNonbondedForce::setComputedValueParameters(int index, const string& name, const string& expression) {
    ASSERT_VALID_INDEX(index, computedValues);
    computedValues[index].name = name;
    computedValues[index].expression = expression;
}

302
int CustomNonbondedForce::addInteractionGroup(const std::set<int>& set1, const std::set<int>& set2) {
303
    for (set<int>::iterator it = set1.begin(); it != set1.end(); ++it)
304
        ASSERT(*it >= 0);
305
    for (set<int>::iterator it = set2.begin(); it != set2.end(); ++it)
306
        ASSERT(*it >= 0);
307
308
309
310
311
312
313
314
315
316
    interactionGroups.push_back(InteractionGroupInfo(set1, set2));
    return interactionGroups.size()-1;
}

void CustomNonbondedForce::getInteractionGroupParameters(int index, std::set<int>& set1, std::set<int>& set2) const {
    ASSERT_VALID_INDEX(index, interactionGroups);
    set1 = interactionGroups[index].set1;
    set2 = interactionGroups[index].set2;
}

peastman's avatar
peastman committed
317
318
void CustomNonbondedForce::setInteractionGroupParameters(int index, const std::set<int>& set1, const std::set<int>& set2) {
    ASSERT_VALID_INDEX(index, interactionGroups);
319
320
321
322
    for (set<int>::iterator it = set1.begin(); it != set1.end(); ++it)
        ASSERT_VALID_INDEX(*it, particles);
    for (set<int>::iterator it = set2.begin(); it != set2.end(); ++it)
        ASSERT_VALID_INDEX(*it, particles);
peastman's avatar
peastman committed
323
324
325
326
    interactionGroups[index].set1 = set1;
    interactionGroups[index].set2 = set2;
}

327
ForceImpl* CustomNonbondedForce::createImpl() const {
328
329
330
331
332
333
    if (numContexts == 0) {
        // Begin tracking changes to particles.
        firstChangedParticle = particles.size();
        lastChangedParticle = -1;
    }
    numContexts++;
334
335
    return new CustomNonbondedForceImpl(*this);
}
336
337

void CustomNonbondedForce::updateParametersInContext(Context& context) {
338
339
340
341
342
343
344
    dynamic_cast<CustomNonbondedForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context), firstChangedParticle, lastChangedParticle);
    if (numContexts == 1) {
        // We just updated the only existing context for this force, so we can reset
        // the tracking of changed particles.
        firstChangedParticle = particles.size();
        lastChangedParticle = -1;
    }
345
}
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381

int CustomNonbondedForce::getNumParticles() const {
    return particles.size();
}

int CustomNonbondedForce::getNumExclusions() const {
    return exclusions.size();
}

int CustomNonbondedForce::getNumPerParticleParameters() const {
    return parameters.size();
}

int CustomNonbondedForce::getNumGlobalParameters() const {
    return globalParameters.size();
}

int CustomNonbondedForce::getNumTabulatedFunctions() const {
    return functions.size();
}

int CustomNonbondedForce::getNumFunctions() const {
    return functions.size();
}

int CustomNonbondedForce::getNumComputedValues() const {
    return computedValues.size();
}

int CustomNonbondedForce::getNumInteractionGroups() const {
    return interactionGroups.size();
}

int CustomNonbondedForce::getNumEnergyParameterDerivatives() const {
    return energyParameterDerivatives.size();
}