ContextImpl.cpp 21.3 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-2016 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
96
    // Validate the list of properties.

    const vector<string>& platformProperties = platform->getPropertyNames();
97
    map<string, string> validatedProperties;
peastman's avatar
peastman committed
98
99
    for (auto& prop : properties) {
        string property = prop.first;
100
101
        if (platform->deprecatedPropertyReplacements.find(property) != platform->deprecatedPropertyReplacements.end())
            property = platform->deprecatedPropertyReplacements[property];
102
        bool valid = false;
peastman's avatar
peastman committed
103
104
        for (auto& p : platformProperties)
            if (p == property) {
105
106
107
108
                valid = true;
                break;
            }
        if (!valid)
peastman's avatar
peastman committed
109
110
            throw OpenMMException("Illegal property name: "+prop.first);
        validatedProperties[property] = prop.second;
111
112
    }
    
113
114
    // Find the list of kernels required.
    
115
    vector<string> kernelNames;
116
    kernelNames.push_back(CalcForcesAndEnergyKernel::Name());
117
    kernelNames.push_back(UpdateStateDataKernel::Name());
118
119
    kernelNames.push_back(ApplyConstraintsKernel::Name());
    kernelNames.push_back(VirtualSitesKernel::Name());
120
    for (int i = 0; i < system.getNumForces(); ++i) {
121
        forceImpls.push_back(system.getForce(i).createImpl());
122
123
124
        vector<string> forceKernels = forceImpls[forceImpls.size()-1]->getKernelNames();
        kernelNames.insert(kernelNames.begin(), forceKernels.begin(), forceKernels.end());
    }
125
    hasInitializedForces = true;
126
127
    vector<string> integratorKernels = integrator.getKernelNames();
    kernelNames.insert(kernelNames.begin(), integratorKernels.begin(), integratorKernels.end());
128
129
130
131
    
    // Select a platform to use.
    
    vector<pair<double, Platform*> > candidatePlatforms;
132
133
134
135
136
    if (platform == NULL) {
        char* defaultPlatform = getenv("OPENMM_DEFAULT_PLATFORM");
        if (defaultPlatform != NULL)
            platform = &Platform::getPlatformByName(string(defaultPlatform));
    }
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
    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;
155
156
157
158
            if (originalContext == NULL)
                platform->contextCreated(*this, validatedProperties);
            else
                platform->linkedContextCreated(*this, *originalContext);
159
160
161
162
163
164
165
166
            break;
        }
        catch (...) {
            if (i > 0)
                continue;
            throw;
        }
    }
167
168
169
}

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

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

211
double ContextImpl::getTime() const {
212
    return updateStateDataKernel.getAs<const UpdateStateDataKernel>().getTime(*this);
213
214
}

215
void ContextImpl::setTime(double t) {
216
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setTime(*this, t);
217
218
}

219
220
221
222
223
224
225
226
long long ContextImpl::getStepCount() const {
    return updateStateDataKernel.getAs<const UpdateStateDataKernel>().getStepCount(*this);
}

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

227
void ContextImpl::getPositions(std::vector<Vec3>& positions) {
228
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPositions(*this, positions);
229
230
231
}

void ContextImpl::setPositions(const std::vector<Vec3>& positions) {
232
    hasSetPositions = true;
233
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPositions(*this, positions);
234
    integrator.stateChanged(State::Positions);
235
236
237
}

void ContextImpl::getVelocities(std::vector<Vec3>& velocities) {
238
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getVelocities(*this, velocities);
239
240
241
}

void ContextImpl::setVelocities(const std::vector<Vec3>& velocities) {
242
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, velocities);
243
    integrator.stateChanged(State::Velocities);
244
245
246
}

void ContextImpl::getForces(std::vector<Vec3>& forces) {
247
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getForces(*this, forces);
248
249
}

250
251
252
253
const std::map<std::string, double>& ContextImpl::getParameters() const {
    return parameters;
}

254
double ContextImpl::getParameter(std::string name) {
255
    if (parameters.find(name) == parameters.end())
256
        throw OpenMMException("Called getParameter() with invalid parameter name: "+name);
257
258
259
    return parameters[name];
}

260
void ContextImpl::setParameter(std::string name, double value) {
261
    if (parameters.find(name) == parameters.end())
262
        throw OpenMMException("Called setParameter() with invalid parameter name: "+name);
263
    parameters[name] = value;
264
    integrator.stateChanged(State::Parameters);
265
266
}

267
268
269
270
void ContextImpl::getEnergyParameterDerivatives(std::map<std::string, double>& derivs) {
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getEnergyParameterDerivatives(*this, derivs);
}

271
void ContextImpl::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) {
272
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPeriodicBoxVectors(*this, a, b, c);
273
274
275
276
277
}

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.");
278
279
280
281
    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.");
282
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, a, b, c);
283
284
}

285
void ContextImpl::applyConstraints(double tol) {
286
287
    if (!hasSetPositions)
        throw OpenMMException("Particle positions have not been set");
288
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().apply(*this, tol);
289
290
}

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

297
void ContextImpl::computeVirtualSites() {
298
    virtualSitesKernel.getAs<VirtualSitesKernel>().computePositions(*this);
299
300
}

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

Peter Eastman's avatar
Peter Eastman committed
318
int& ContextImpl::getLastForceGroups() {
319
320
321
    return lastForceGroups;
}

322
double ContextImpl::calcKineticEnergy() {
323
    return integrator.computeKineticEnergy();
324
325
}

326
327
bool ContextImpl::updateContextState() {
    bool forcesInvalid = false;
peastman's avatar
peastman committed
328
    for (auto force : forceImpls)
329
330
        force->updateContextState(*this, forcesInvalid);
    return forcesInvalid;
331
332
}

333
334
335
336
const vector<ForceImpl*>& ContextImpl::getForceImpls() const {
    return forceImpls;
}

337
338
339
340
vector<ForceImpl*>& ContextImpl::getForceImpls() {
    return forceImpls;
}

341
void* ContextImpl::getPlatformData() {
342
343
344
    return platformData;
}

345
346
347
348
const void* ContextImpl::getPlatformData() const {
    return platformData;
}

349
void ContextImpl::setPlatformData(void* data) {
350
351
    platformData = data;
}
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367

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
368
369
    for (auto force : forceImpls) {
        vector<pair<int, int> > forceBonds = force->getBondedParticles();
370
371
        bonds.insert(bonds.end(), forceBonds.begin(), forceBonds.end());
    }
372
373
374
375
376
377
378
    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)));
        }
    }
379
380
381
382
383

    // 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
384
385
386
    for (auto& bond : bonds) {
        particleBonds[bond.first].push_back(bond.second);
        particleBonds[bond.second].push_back(bond.first);
387
388
    }

389
    // Now identify particles by which molecule they belong to.
390

391
392
393
394
395
396
397
398
399
    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.
    
400
401
402
    vector<int> particleMolecule(numParticles, -1);
    int numMolecules = 0;
    for (int i = 0; i < numParticles; i++)
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
        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);
434
435
436
437
438
    for (int i = 0; i < numParticles; i++)
        molecules[particleMolecule[i]].push_back(i);
    return molecules;
}

Peter Eastman's avatar
Peter Eastman committed
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
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
454
    stream.write(CHECKPOINT_MAGIC_BYTES, sizeof(CHECKPOINT_MAGIC_BYTES)/sizeof(CHECKPOINT_MAGIC_BYTES[0]));
Peter Eastman's avatar
Peter Eastman committed
455
456
457
458
459
    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
460
461
462
    for (auto& param : parameters) {
        writeString(stream, param.first);
        stream.write((char*) &param.second, sizeof(double));
Peter Eastman's avatar
Peter Eastman committed
463
    }
464
    updateStateDataKernel.getAs<UpdateStateDataKernel>().createCheckpoint(*this, stream);
465
    integrator.createCheckpoint(stream);
Peter Eastman's avatar
Peter Eastman committed
466
467
468
469
    stream.flush();
}

void ContextImpl::loadCheckpoint(istream& stream) {
Robert McGibbon's avatar
Robert McGibbon committed
470
471
472
473
    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
474
475
        throw OpenMMException("loadCheckpoint: Checkpoint header was not correct");

Peter Eastman's avatar
Peter Eastman committed
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
    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;
    }
491
    updateStateDataKernel.getAs<UpdateStateDataKernel>().loadCheckpoint(*this, stream);
492
    integrator.loadCheckpoint(stream);
493
    hasSetPositions = true;
494
495
496
497
    integrator.stateChanged(State::Positions);
    integrator.stateChanged(State::Velocities);
    integrator.stateChanged(State::Parameters);
    integrator.stateChanged(State::Energy);
Peter Eastman's avatar
Peter Eastman committed
498
}
499
500
501
502

void ContextImpl::systemChanged() {
    integrator.stateChanged(State::Energy);
}
503
504
505
506

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