ContextImpl.cpp 21.7 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-2022 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
 * 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.                                     *
 * -------------------------------------------------------------------------- */

32
33
34
35
36
37
#include "openmm/Force.h"
#include "openmm/Integrator.h"
#include "openmm/OpenMMException.h"
#include "openmm/System.h"
#include "openmm/kernels.h"
#include "openmm/internal/ForceImpl.h"
38
#include "openmm/internal/ContextImpl.h"
39
#include "openmm/State.h"
40
#include "openmm/VirtualSite.h"
Peter Eastman's avatar
Peter Eastman committed
41
#include "openmm/Context.h"
42
#include <algorithm>
43
#include <cmath>
44
#include <cstdlib>
Peter Eastman's avatar
Peter Eastman committed
45
#include <iostream>
46
#include <map>
47
#include <utility>
48
#include <vector>
Robert McGibbon's avatar
Robert McGibbon committed
49
#include <string.h>
50
51

using namespace OpenMM;
Peter Eastman's avatar
Peter Eastman committed
52
using namespace std;
Robert McGibbon's avatar
Robert McGibbon committed
53
const static char CHECKPOINT_MAGIC_BYTES[] = "OpenMM Binary Checkpoint\n";
Robert McGibbon's avatar
Robert McGibbon committed
54

55

56
ContextImpl::ContextImpl(Context& owner, const System& system, Integrator& integrator, Platform* platform, const map<string, string>& properties, ContextImpl* originalContext) :
57
58
        owner(owner), system(system), integrator(integrator), hasInitializedForces(false), hasSetPositions(false), integratorIsDeleted(false),
        lastForceGroups(-1), platform(platform), platformData(NULL) {
59
60
    int numParticles = system.getNumParticles();
    if (numParticles == 0)
61
        throw OpenMMException("Cannot create a Context for a System with no particles");
62
63
64
    
    // Check for errors in virtual sites and massless particles.
    
65
    for (int i = 0; i < numParticles; i++) {
66
67
68
69
70
71
72
73
74
        if (system.isVirtualSite(i)) {
            if (system.getParticleMass(i) != 0.0)
                throw OpenMMException("Virtual site has nonzero mass");
            const VirtualSite& site = system.getVirtualSite(i);
            for (int j = 0; j < site.getNumParticles(); j++)
                if (system.isVirtualSite(site.getParticle(j)))
                    throw OpenMMException("A virtual site cannot depend on another virtual site");
        }
    }
75
    set<pair<int, int> > constraintAtoms;
76
77
78
79
    for (int i = 0; i < system.getNumConstraints(); i++) {
        int particle1, particle2;
        double distance;
        system.getConstraintParameters(i, particle1, particle2, distance);
80
81
82
83
        if (particle1 == particle2)
            throw OpenMMException("A constraint cannot connect a particle to itself");
        if (particle1 < 0 || particle2 < 0 || particle1 >= numParticles || particle2 >= numParticles)
            throw OpenMMException("Illegal particle index in constraint");
84
85
86
        double mass1 = system.getParticleMass(particle1);
        double mass2 = system.getParticleMass(particle2);
        if ((mass1 == 0.0 && mass2 != 0.0) || (mass2 == 0.0 && mass1 != 0.0))
87
            throw OpenMMException("A constraint cannot involve a massless particle");
88
89
90
91
        pair<int, int> atoms = make_pair(min(particle1, particle2), max(particle1, particle2));
        if (constraintAtoms.find(atoms) != constraintAtoms.end())
            throw OpenMMException("The System has two constraints between the same atoms.  This will produce a singular constraint matrix.");
        constraintAtoms.insert(atoms);
92
93
    }
    
94
95
    // Validate the list of properties.

96
    map<string, string> validatedProperties;
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
    if (platform != NULL) {
        const vector<string>& platformProperties = platform->getPropertyNames();
        for (auto& prop : properties) {
            string property = prop.first;
            if (platform->deprecatedPropertyReplacements.find(property) != platform->deprecatedPropertyReplacements.end())
                property = platform->deprecatedPropertyReplacements[property];
            bool valid = false;
            for (auto& p : platformProperties)
                if (p == property) {
                    valid = true;
                    break;
                }
            if (!valid)
                throw OpenMMException("Illegal property name: "+prop.first);
            validatedProperties[property] = prop.second;
        }
113
    }
114
115
116
117
    else {
        // There can't be any platform-specific properties if there's no platform
    }

118
119
    // Find the list of kernels required.
    
120
    vector<string> kernelNames;
121
    kernelNames.push_back(CalcForcesAndEnergyKernel::Name());
122
    kernelNames.push_back(UpdateStateDataKernel::Name());
123
124
    kernelNames.push_back(ApplyConstraintsKernel::Name());
    kernelNames.push_back(VirtualSitesKernel::Name());
125
    for (int i = 0; i < system.getNumForces(); ++i) {
126
        forceImpls.push_back(system.getForce(i).createImpl());
127
128
129
        vector<string> forceKernels = forceImpls[forceImpls.size()-1]->getKernelNames();
        kernelNames.insert(kernelNames.begin(), forceKernels.begin(), forceKernels.end());
    }
130
    hasInitializedForces = true;
131
132
    vector<string> integratorKernels = integrator.getKernelNames();
    kernelNames.insert(kernelNames.begin(), integratorKernels.begin(), integratorKernels.end());
133
134
135
136
    
    // Select a platform to use.
    
    vector<pair<double, Platform*> > candidatePlatforms;
137
138
139
140
141
    if (platform == NULL) {
        char* defaultPlatform = getenv("OPENMM_DEFAULT_PLATFORM");
        if (defaultPlatform != NULL)
            platform = &Platform::getPlatformByName(string(defaultPlatform));
    }
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
    if (platform == NULL) {
        for (int i = 0; i < Platform::getNumPlatforms(); i++) {
            Platform& p = Platform::getPlatform(i);
            if (p.supportsKernels(kernelNames))
                candidatePlatforms.push_back(make_pair(p.getSpeed(), &p));
        }
        if (candidatePlatforms.size() == 0)
            throw OpenMMException("No Platform supports all the requested kernels");
        sort(candidatePlatforms.begin(), candidatePlatforms.end());
    }
    else {
        if (!platform->supportsKernels(kernelNames))
            throw OpenMMException("Specified a Platform for a Context which does not support all required kernels");
        candidatePlatforms.push_back(make_pair(platform->getSpeed(), platform));
    }
    for (int i = candidatePlatforms.size()-1; i >= 0; i--) {
        try {
            this->platform = platform = candidatePlatforms[i].second;
160
161
162
163
            if (originalContext == NULL)
                platform->contextCreated(*this, validatedProperties);
            else
                platform->linkedContextCreated(*this, *originalContext);
164
165
166
167
168
169
170
171
            break;
        }
        catch (...) {
            if (i > 0)
                continue;
            throw;
        }
    }
172
173
174
}

void ContextImpl::initialize() {
175
176
    // Create and initialize kernels and other objects.
    
177
    initializeForcesKernel = platform->createKernel(CalcForcesAndEnergyKernel::Name(), *this);
178
    initializeForcesKernel.getAs<CalcForcesAndEnergyKernel>().initialize(system);
179
    updateStateDataKernel = platform->createKernel(UpdateStateDataKernel::Name(), *this);
180
    updateStateDataKernel.getAs<UpdateStateDataKernel>().initialize(system);
181
    applyConstraintsKernel = platform->createKernel(ApplyConstraintsKernel::Name(), *this);
182
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().initialize(system);
183
    virtualSitesKernel = platform->createKernel(VirtualSitesKernel::Name(), *this);
184
    virtualSitesKernel.getAs<VirtualSitesKernel>().initialize(system);
185
186
    Vec3 periodicBoxVectors[3];
    system.getDefaultPeriodicBoxVectors(periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
187
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
188
    for (size_t i = 0; i < forceImpls.size(); ++i) {
189
        forceImpls[i]->initialize(*this);
190
191
192
        map<string, double> forceParameters = forceImpls[i]->getDefaultParameters();
        parameters.insert(forceParameters.begin(), forceParameters.end());
    }
193
    integrator.initialize(*this);
194
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, vector<Vec3>(system.getNumParticles()));
195
196
}

197
ContextImpl::~ContextImpl() {
peastman's avatar
peastman committed
198
199
    for (auto force : forceImpls)
        delete force;
200
201
202
203
204
205
206
    
    // Make sure all kernels get properly deleted before contextDestroyed() is called.
    
    initializeForcesKernel = Kernel();
    updateStateDataKernel = Kernel();
    applyConstraintsKernel = Kernel();
    virtualSitesKernel = Kernel();
207
208
209
210
211
212
    if (!integratorIsDeleted) {
        // The Context is being deleted before the Integrator, so call cleanup() on it now.
        
        integrator.cleanup();
        integrator.context = NULL;
    }
213
    platform->contextDestroyed(*this);
214
215
}

216
double ContextImpl::getTime() const {
217
    return updateStateDataKernel.getAs<const UpdateStateDataKernel>().getTime(*this);
218
219
}

220
void ContextImpl::setTime(double t) {
221
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setTime(*this, t);
222
223
}

224
225
226
227
228
229
230
231
long long ContextImpl::getStepCount() const {
    return updateStateDataKernel.getAs<const UpdateStateDataKernel>().getStepCount(*this);
}

void ContextImpl::setStepCount(long long count) {
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setStepCount(*this, count);
}

232
void ContextImpl::getPositions(std::vector<Vec3>& positions) {
233
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPositions(*this, positions);
234
235
236
}

void ContextImpl::setPositions(const std::vector<Vec3>& positions) {
237
    hasSetPositions = true;
238
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPositions(*this, positions);
239
    integrator.stateChanged(State::Positions);
240
241
242
}

void ContextImpl::getVelocities(std::vector<Vec3>& velocities) {
243
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getVelocities(*this, velocities);
244
245
246
}

void ContextImpl::setVelocities(const std::vector<Vec3>& velocities) {
247
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, velocities);
248
    integrator.stateChanged(State::Velocities);
249
250
251
}

void ContextImpl::getForces(std::vector<Vec3>& forces) {
252
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getForces(*this, forces);
253
254
}

255
256
257
258
const std::map<std::string, double>& ContextImpl::getParameters() const {
    return parameters;
}

259
double ContextImpl::getParameter(std::string name) {
260
    if (parameters.find(name) == parameters.end())
261
        throw OpenMMException("Called getParameter() with invalid parameter name: "+name);
262
263
264
    return parameters[name];
}

265
void ContextImpl::setParameter(std::string name, double value) {
266
    if (parameters.find(name) == parameters.end())
267
        throw OpenMMException("Called setParameter() with invalid parameter name: "+name);
268
    parameters[name] = value;
269
    integrator.stateChanged(State::Parameters);
270
271
}

272
273
274
275
void ContextImpl::getEnergyParameterDerivatives(std::map<std::string, double>& derivs) {
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getEnergyParameterDerivatives(*this, derivs);
}

276
void ContextImpl::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) {
277
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPeriodicBoxVectors(*this, a, b, c);
278
279
280
281
282
}

void ContextImpl::setPeriodicBoxVectors(const Vec3& a, const Vec3& b, const Vec3& c) {
    if (a[1] != 0.0 || a[2] != 0.0)
        throw OpenMMException("First periodic box vector must be parallel to x.");
283
284
285
286
    if (b[2] != 0.0)
        throw OpenMMException("Second periodic box vector must be in the x-y plane.");
    if (a[0] <= 0.0 || b[1] <= 0.0 || c[2] <= 0.0 || a[0] < 2*fabs(b[0]) || a[0] < 2*fabs(c[0]) || b[1] < 2*fabs(c[1]))
        throw OpenMMException("Periodic box vectors must be in reduced form.");
287
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, a, b, c);
288
289
}

290
void ContextImpl::applyConstraints(double tol) {
291
292
    if (!hasSetPositions)
        throw OpenMMException("Particle positions have not been set");
293
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().apply(*this, tol);
294
295
}

296
void ContextImpl::applyVelocityConstraints(double tol) {
297
298
    if (!hasSetPositions)
        throw OpenMMException("Particle positions have not been set");
299
300
301
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().applyToVelocities(*this, tol);
}

302
void ContextImpl::computeVirtualSites() {
303
    virtualSitesKernel.getAs<VirtualSitesKernel>().computePositions(*this);
304
305
}

306
double ContextImpl::calcForcesAndEnergy(bool includeForces, bool includeEnergy, int groups) {
307
308
    if (!hasSetPositions)
        throw OpenMMException("Particle positions have not been set");
309
    lastForceGroups = groups;
310
    CalcForcesAndEnergyKernel& kernel = initializeForcesKernel.getAs<CalcForcesAndEnergyKernel>();
311
312
313
    while (true) {
        double energy = 0.0;
        kernel.beginComputation(*this, includeForces, includeEnergy, groups);
peastman's avatar
peastman committed
314
315
        for (auto force : forceImpls)
            energy += force->calcForcesAndEnergy(*this, includeForces, includeEnergy, groups);
316
317
318
319
320
        bool valid = true;
        energy += kernel.finishComputation(*this, includeForces, includeEnergy, groups, valid);
        if (valid)
            return energy;
    }
321
322
}

Peter Eastman's avatar
Peter Eastman committed
323
int& ContextImpl::getLastForceGroups() {
324
325
326
    return lastForceGroups;
}

327
double ContextImpl::calcKineticEnergy() {
328
    return integrator.computeKineticEnergy();
329
330
}

331
332
333
334
void ContextImpl::computeShiftedVelocities(double timeShift, std::vector<Vec3>& velocities) {
    updateStateDataKernel.getAs<UpdateStateDataKernel>().computeShiftedVelocities(*this, timeShift, velocities);
}

335
336
bool ContextImpl::updateContextState() {
    bool forcesInvalid = false;
peastman's avatar
peastman committed
337
    for (auto force : forceImpls)
338
339
        force->updateContextState(*this, forcesInvalid);
    return forcesInvalid;
340
341
}

342
343
344
345
const vector<ForceImpl*>& ContextImpl::getForceImpls() const {
    return forceImpls;
}

346
347
348
349
vector<ForceImpl*>& ContextImpl::getForceImpls() {
    return forceImpls;
}

350
void* ContextImpl::getPlatformData() {
351
352
353
    return platformData;
}

354
355
356
357
const void* ContextImpl::getPlatformData() const {
    return platformData;
}

358
void ContextImpl::setPlatformData(void* data) {
359
360
    platformData = data;
}
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376

const vector<vector<int> >& ContextImpl::getMolecules() const {
    if (!hasInitializedForces)
        throw OpenMMException("ContextImpl: getMolecules() cannot be called until all ForceImpls have been initialized");
    if (molecules.size() > 0 || system.getNumParticles() == 0)
        return molecules;

    // First make a list of bonds and constraints.

    vector<pair<int, int> > bonds;
    for (int i = 0; i < system.getNumConstraints(); i++) {
        int particle1, particle2;
        double distance;
        system.getConstraintParameters(i, particle1, particle2, distance);
        bonds.push_back(std::make_pair(particle1, particle2));
    }
peastman's avatar
peastman committed
377
378
    for (auto force : forceImpls) {
        vector<pair<int, int> > forceBonds = force->getBondedParticles();
379
380
        bonds.insert(bonds.end(), forceBonds.begin(), forceBonds.end());
    }
381
382
383
384
385
386
387
    for (int i = 0; i < system.getNumParticles(); i++) {
        if (system.isVirtualSite(i)) {
            const VirtualSite& site = system.getVirtualSite(i);
            for (int j = 0; j < site.getNumParticles(); j++)
                bonds.push_back(std::make_pair(i, site.getParticle(j)));
        }
    }
388
389
390
391
392

    // Make a list of every other particle to which each particle is connected

    int numParticles = system.getNumParticles();
    vector<vector<int> > particleBonds(numParticles);
peastman's avatar
peastman committed
393
394
395
    for (auto& bond : bonds) {
        particleBonds[bond.first].push_back(bond.second);
        particleBonds[bond.second].push_back(bond.first);
396
397
    }

398
    // Now identify particles by which molecule they belong to.
399

400
401
402
403
404
405
406
407
408
    molecules = findMolecules(numParticles, particleBonds);
    return molecules;
}

vector<vector<int> > ContextImpl::findMolecules(int numParticles, vector<vector<int> >& particleBonds) {
    // This is essentially a recursive algorithm, but it is reformulated as a loop to avoid
    // stack overflows.  It selects a particle, marks it as a new molecule, then recursively
    // marks every particle bonded to it as also being in that molecule.
    
409
410
411
    vector<int> particleMolecule(numParticles, -1);
    int numMolecules = 0;
    for (int i = 0; i < numParticles; i++)
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
        if (particleMolecule[i] == -1) {
            // Start a new molecule.
            
            vector<int> particleStack;
            vector<int> neighborStack;
            particleStack.push_back(i);
            neighborStack.push_back(0);
            int molecule = numMolecules++;
            
            // Recursively tag all the bonded particles.
            
            while (particleStack.size() > 0) {
                int particle = particleStack.back();
                particleMolecule[particle] = molecule;
                int& neighbor = neighborStack.back();
                while (neighbor < particleBonds[particle].size() && particleMolecule[particleBonds[particle][neighbor]] != -1)
                    neighbor++;
                if (neighbor < particleBonds[particle].size()) {
                    particleStack.push_back(particleBonds[particle][neighbor]);
                    neighborStack.push_back(0);
                }
                else {
                    particleStack.pop_back();
                    neighborStack.pop_back();
                }
            }
        }
    
    // Build the final output vector.
    
    vector<vector<int> > molecules(numMolecules);
443
444
445
446
447
    for (int i = 0; i < numParticles; i++)
        molecules[particleMolecule[i]].push_back(i);
    return molecules;
}

Peter Eastman's avatar
Peter Eastman committed
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
static void writeString(ostream& stream, string str) {
    int length = str.size();
    stream.write((char*) &length, sizeof(int));
    stream.write((char*) &str[0], length);
}

static string readString(istream& stream) {
    int length;
    stream.read((char*) &length, sizeof(int));
    string str(length, ' ');
    stream.read((char*) &str[0], length);
    return str;
}

void ContextImpl::createCheckpoint(ostream& stream) {
Robert McGibbon's avatar
Robert McGibbon committed
463
    stream.write(CHECKPOINT_MAGIC_BYTES, sizeof(CHECKPOINT_MAGIC_BYTES)/sizeof(CHECKPOINT_MAGIC_BYTES[0]));
Peter Eastman's avatar
Peter Eastman committed
464
465
466
467
468
    writeString(stream, getPlatform().getName());
    int numParticles = getSystem().getNumParticles();
    stream.write((char*) &numParticles, sizeof(int));
    int numParameters = parameters.size();
    stream.write((char*) &numParameters, sizeof(int));
peastman's avatar
peastman committed
469
470
471
    for (auto& param : parameters) {
        writeString(stream, param.first);
        stream.write((char*) &param.second, sizeof(double));
Peter Eastman's avatar
Peter Eastman committed
472
    }
473
    updateStateDataKernel.getAs<UpdateStateDataKernel>().createCheckpoint(*this, stream);
474
    integrator.createCheckpoint(stream);
Peter Eastman's avatar
Peter Eastman committed
475
476
477
478
    stream.flush();
}

void ContextImpl::loadCheckpoint(istream& stream) {
Robert McGibbon's avatar
Robert McGibbon committed
479
480
481
482
    static const int magiclength = sizeof(CHECKPOINT_MAGIC_BYTES)/sizeof(CHECKPOINT_MAGIC_BYTES[0]);
    char magicbytes[magiclength];
    stream.read(magicbytes, magiclength);
    if (memcmp(magicbytes, CHECKPOINT_MAGIC_BYTES, magiclength) != 0)
Robert McGibbon's avatar
Robert McGibbon committed
483
484
        throw OpenMMException("loadCheckpoint: Checkpoint header was not correct");

Peter Eastman's avatar
Peter Eastman committed
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
    string platformName = readString(stream);
    if (platformName != getPlatform().getName())
        throw OpenMMException("loadCheckpoint: Checkpoint was created with a different Platform: "+platformName);
    int numParticles;
    stream.read((char*) &numParticles, sizeof(int));
    if (numParticles != getSystem().getNumParticles())
        throw OpenMMException("loadCheckpoint: Checkpoint contains the wrong number of particles");
    int numParameters;
    stream.read((char*) &numParameters, sizeof(int));
    for (int i = 0; i < numParameters; i++) {
        string name = readString(stream);
        double value;
        stream.read((char*) &value, sizeof(double));
        parameters[name] = value;
    }
500
    updateStateDataKernel.getAs<UpdateStateDataKernel>().loadCheckpoint(*this, stream);
501
    integrator.loadCheckpoint(stream);
502
    hasSetPositions = true;
503
504
505
506
    integrator.stateChanged(State::Positions);
    integrator.stateChanged(State::Velocities);
    integrator.stateChanged(State::Parameters);
    integrator.stateChanged(State::Energy);
Peter Eastman's avatar
Peter Eastman committed
507
}
508
509
510
511

void ContextImpl::systemChanged() {
    integrator.stateChanged(State::Energy);
}
512
513
514
515

Context* ContextImpl::createLinkedContext(const System& system, Integrator& integrator) {
    return new Context(system, integrator, *this);
}