ContextImpl.cpp 20.4 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) :
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
            platform->contextCreated(*this, validatedProperties);
156
157
158
159
160
161
162
163
            break;
        }
        catch (...) {
            if (i > 0)
                continue;
            throw;
        }
    }
164
165
166
    
    // Create and initialize kernels and other objects.
    
167
    initializeForcesKernel = platform->createKernel(CalcForcesAndEnergyKernel::Name(), *this);
168
    initializeForcesKernel.getAs<CalcForcesAndEnergyKernel>().initialize(system);
169
    updateStateDataKernel = platform->createKernel(UpdateStateDataKernel::Name(), *this);
170
    updateStateDataKernel.getAs<UpdateStateDataKernel>().initialize(system);
171
    applyConstraintsKernel = platform->createKernel(ApplyConstraintsKernel::Name(), *this);
172
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().initialize(system);
173
    virtualSitesKernel = platform->createKernel(VirtualSitesKernel::Name(), *this);
174
    virtualSitesKernel.getAs<VirtualSitesKernel>().initialize(system);
175
176
    Vec3 periodicBoxVectors[3];
    system.getDefaultPeriodicBoxVectors(periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
177
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
178
    for (size_t i = 0; i < forceImpls.size(); ++i) {
179
        forceImpls[i]->initialize(*this);
180
181
182
        map<string, double> forceParameters = forceImpls[i]->getDefaultParameters();
        parameters.insert(forceParameters.begin(), forceParameters.end());
    }
183
    integrator.initialize(*this);
184
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, vector<Vec3>(numParticles));
185
186
}

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

206
double ContextImpl::getTime() const {
207
    return updateStateDataKernel.getAs<const UpdateStateDataKernel>().getTime(*this);
208
209
}

210
void ContextImpl::setTime(double t) {
211
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setTime(*this, t);
212
213
214
}

void ContextImpl::getPositions(std::vector<Vec3>& positions) {
215
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPositions(*this, positions);
216
217
218
}

void ContextImpl::setPositions(const std::vector<Vec3>& positions) {
219
    hasSetPositions = true;
220
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPositions(*this, positions);
221
    integrator.stateChanged(State::Positions);
222
223
224
}

void ContextImpl::getVelocities(std::vector<Vec3>& velocities) {
225
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getVelocities(*this, velocities);
226
227
228
}

void ContextImpl::setVelocities(const std::vector<Vec3>& velocities) {
229
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, velocities);
230
    integrator.stateChanged(State::Velocities);
231
232
233
}

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

237
238
239
240
const std::map<std::string, double>& ContextImpl::getParameters() const {
    return parameters;
}

241
double ContextImpl::getParameter(std::string name) {
242
    if (parameters.find(name) == parameters.end())
243
        throw OpenMMException("Called getParameter() with invalid parameter name: "+name);
244
245
246
    return parameters[name];
}

247
void ContextImpl::setParameter(std::string name, double value) {
248
    if (parameters.find(name) == parameters.end())
249
        throw OpenMMException("Called setParameter() with invalid parameter name: "+name);
250
    parameters[name] = value;
251
    integrator.stateChanged(State::Parameters);
252
253
}

254
255
256
257
void ContextImpl::getEnergyParameterDerivatives(std::map<std::string, double>& derivs) {
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getEnergyParameterDerivatives(*this, derivs);
}

258
void ContextImpl::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) {
259
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPeriodicBoxVectors(*this, a, b, c);
260
261
262
263
264
}

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.");
265
266
267
268
    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.");
269
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, a, b, c);
270
271
}

272
void ContextImpl::applyConstraints(double tol) {
273
274
    if (!hasSetPositions)
        throw OpenMMException("Particle positions have not been set");
275
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().apply(*this, tol);
276
277
}

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

284
void ContextImpl::computeVirtualSites() {
285
    virtualSitesKernel.getAs<VirtualSitesKernel>().computePositions(*this);
286
287
}

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

Peter Eastman's avatar
Peter Eastman committed
305
int& ContextImpl::getLastForceGroups() {
306
307
308
    return lastForceGroups;
}

309
double ContextImpl::calcKineticEnergy() {
310
    return integrator.computeKineticEnergy();
311
312
}

313
void ContextImpl::updateContextState() {
peastman's avatar
peastman committed
314
315
    for (auto force : forceImpls)
        force->updateContextState(*this);
316
317
}

318
319
320
321
const vector<ForceImpl*>& ContextImpl::getForceImpls() const {
    return forceImpls;
}

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

326
void* ContextImpl::getPlatformData() {
327
328
329
    return platformData;
}

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

334
void ContextImpl::setPlatformData(void* data) {
335
336
    platformData = data;
}
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352

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

    // 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
369
370
371
    for (auto& bond : bonds) {
        particleBonds[bond.first].push_back(bond.second);
        particleBonds[bond.second].push_back(bond.first);
372
373
    }

374
    // Now identify particles by which molecule they belong to.
375

376
377
378
379
380
381
382
383
384
    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.
    
385
386
387
    vector<int> particleMolecule(numParticles, -1);
    int numMolecules = 0;
    for (int i = 0; i < numParticles; i++)
388
389
390
391
392
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
        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);
419
420
421
422
423
    for (int i = 0; i < numParticles; i++)
        molecules[particleMolecule[i]].push_back(i);
    return molecules;
}

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

void ContextImpl::loadCheckpoint(istream& stream) {
Robert McGibbon's avatar
Robert McGibbon committed
454
455
456
457
    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
458
459
        throw OpenMMException("loadCheckpoint: Checkpoint header was not correct");

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

void ContextImpl::systemChanged() {
    integrator.stateChanged(State::Energy);
}