ContextImpl.cpp 20.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-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
}

void ContextImpl::getPositions(std::vector<Vec3>& positions) {
220
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPositions(*this, positions);
221
222
223
}

void ContextImpl::setPositions(const std::vector<Vec3>& positions) {
224
    hasSetPositions = true;
225
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPositions(*this, positions);
226
    integrator.stateChanged(State::Positions);
227
228
229
}

void ContextImpl::getVelocities(std::vector<Vec3>& velocities) {
230
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getVelocities(*this, velocities);
231
232
233
}

void ContextImpl::setVelocities(const std::vector<Vec3>& velocities) {
234
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, velocities);
235
    integrator.stateChanged(State::Velocities);
236
237
238
}

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

242
243
244
245
const std::map<std::string, double>& ContextImpl::getParameters() const {
    return parameters;
}

246
double ContextImpl::getParameter(std::string name) {
247
    if (parameters.find(name) == parameters.end())
248
        throw OpenMMException("Called getParameter() with invalid parameter name: "+name);
249
250
251
    return parameters[name];
}

252
void ContextImpl::setParameter(std::string name, double value) {
253
    if (parameters.find(name) == parameters.end())
254
        throw OpenMMException("Called setParameter() with invalid parameter name: "+name);
255
    parameters[name] = value;
256
    integrator.stateChanged(State::Parameters);
257
258
}

259
260
261
262
void ContextImpl::getEnergyParameterDerivatives(std::map<std::string, double>& derivs) {
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getEnergyParameterDerivatives(*this, derivs);
}

263
void ContextImpl::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) {
264
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPeriodicBoxVectors(*this, a, b, c);
265
266
267
268
269
}

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.");
270
271
272
273
    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.");
274
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, a, b, c);
275
276
}

277
void ContextImpl::applyConstraints(double tol) {
278
279
    if (!hasSetPositions)
        throw OpenMMException("Particle positions have not been set");
280
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().apply(*this, tol);
281
282
}

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

289
void ContextImpl::computeVirtualSites() {
290
    virtualSitesKernel.getAs<VirtualSitesKernel>().computePositions(*this);
291
292
}

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

Peter Eastman's avatar
Peter Eastman committed
310
int& ContextImpl::getLastForceGroups() {
311
312
313
    return lastForceGroups;
}

314
double ContextImpl::calcKineticEnergy() {
315
    return integrator.computeKineticEnergy();
316
317
}

318
void ContextImpl::updateContextState() {
peastman's avatar
peastman committed
319
320
    for (auto force : forceImpls)
        force->updateContextState(*this);
321
322
}

323
324
325
326
const vector<ForceImpl*>& ContextImpl::getForceImpls() const {
    return forceImpls;
}

327
328
329
330
vector<ForceImpl*>& ContextImpl::getForceImpls() {
    return forceImpls;
}

331
void* ContextImpl::getPlatformData() {
332
333
334
    return platformData;
}

335
336
337
338
const void* ContextImpl::getPlatformData() const {
    return platformData;
}

339
void ContextImpl::setPlatformData(void* data) {
340
341
    platformData = data;
}
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357

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
358
359
    for (auto force : forceImpls) {
        vector<pair<int, int> > forceBonds = force->getBondedParticles();
360
361
        bonds.insert(bonds.end(), forceBonds.begin(), forceBonds.end());
    }
362
363
364
365
366
367
368
    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)));
        }
    }
369
370
371
372
373

    // 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
374
375
376
    for (auto& bond : bonds) {
        particleBonds[bond.first].push_back(bond.second);
        particleBonds[bond.second].push_back(bond.first);
377
378
    }

379
    // Now identify particles by which molecule they belong to.
380

381
382
383
384
385
386
387
388
389
    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.
    
390
391
392
    vector<int> particleMolecule(numParticles, -1);
    int numMolecules = 0;
    for (int i = 0; i < numParticles; i++)
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
        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);
424
425
426
427
428
    for (int i = 0; i < numParticles; i++)
        molecules[particleMolecule[i]].push_back(i);
    return molecules;
}

Peter Eastman's avatar
Peter Eastman committed
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
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
444
    stream.write(CHECKPOINT_MAGIC_BYTES, sizeof(CHECKPOINT_MAGIC_BYTES)/sizeof(CHECKPOINT_MAGIC_BYTES[0]));
Peter Eastman's avatar
Peter Eastman committed
445
446
447
448
449
    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
450
451
452
    for (auto& param : parameters) {
        writeString(stream, param.first);
        stream.write((char*) &param.second, sizeof(double));
Peter Eastman's avatar
Peter Eastman committed
453
    }
454
    updateStateDataKernel.getAs<UpdateStateDataKernel>().createCheckpoint(*this, stream);
Peter Eastman's avatar
Peter Eastman committed
455
456
457
458
    stream.flush();
}

void ContextImpl::loadCheckpoint(istream& stream) {
Robert McGibbon's avatar
Robert McGibbon committed
459
460
461
462
    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
463
464
        throw OpenMMException("loadCheckpoint: Checkpoint header was not correct");

Peter Eastman's avatar
Peter Eastman committed
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
    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;
    }
480
    updateStateDataKernel.getAs<UpdateStateDataKernel>().loadCheckpoint(*this, stream);
481
    hasSetPositions = true;
Peter Eastman's avatar
Peter Eastman committed
482
}
483
484
485
486

void ContextImpl::systemChanged() {
    integrator.stateChanged(State::Energy);
}
487
488
489
490

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