ContextImpl.cpp 19.6 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
    if (system.getNumParticles() == 0)
        throw OpenMMException("Cannot create a Context for a System with no particles");
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
    
    // Check for errors in virtual sites and massless particles.
    
    for (int i = 0; i < system.getNumParticles(); i++) {
        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");
        }
    }
    for (int i = 0; i < system.getNumConstraints(); i++) {
        int particle1, particle2;
        double distance;
        system.getConstraintParameters(i, particle1, particle2, distance);
78
79
80
        double mass1 = system.getParticleMass(particle1);
        double mass2 = system.getParticleMass(particle2);
        if ((mass1 == 0.0 && mass2 != 0.0) || (mass2 == 0.0 && mass1 != 0.0))
81
82
83
            throw OpenMMException("A constraint cannot involve a massless particle");
    }
    
84
85
86
    // Validate the list of properties.

    const vector<string>& platformProperties = platform->getPropertyNames();
87
    map<string, string> validatedProperties;
88
    for (map<string, string>::const_iterator iter = properties.begin(); iter != properties.end(); ++iter) {
89
90
91
        string property = iter->first;
        if (platform->deprecatedPropertyReplacements.find(property) != platform->deprecatedPropertyReplacements.end())
            property = platform->deprecatedPropertyReplacements[property];
92
93
        bool valid = false;
        for (int i = 0; i < (int) platformProperties.size(); i++)
94
            if (platformProperties[i] == property) {
95
96
97
98
99
                valid = true;
                break;
            }
        if (!valid)
            throw OpenMMException("Illegal property name: "+iter->first);
100
        validatedProperties[property] = iter->second;
101
102
    }
    
103
104
    // Find the list of kernels required.
    
105
    vector<string> kernelNames;
106
    kernelNames.push_back(CalcForcesAndEnergyKernel::Name());
107
    kernelNames.push_back(UpdateStateDataKernel::Name());
108
109
    kernelNames.push_back(ApplyConstraintsKernel::Name());
    kernelNames.push_back(VirtualSitesKernel::Name());
110
    for (int i = 0; i < system.getNumForces(); ++i) {
111
        forceImpls.push_back(system.getForce(i).createImpl());
112
113
        map<string, double> forceParameters = forceImpls[forceImpls.size()-1]->getDefaultParameters();
        parameters.insert(forceParameters.begin(), forceParameters.end());
114
115
116
        vector<string> forceKernels = forceImpls[forceImpls.size()-1]->getKernelNames();
        kernelNames.insert(kernelNames.begin(), forceKernels.begin(), forceKernels.end());
    }
117
    hasInitializedForces = true;
118
119
    vector<string> integratorKernels = integrator.getKernelNames();
    kernelNames.insert(kernelNames.begin(), integratorKernels.begin(), integratorKernels.end());
120
121
122
123
    
    // Select a platform to use.
    
    vector<pair<double, Platform*> > candidatePlatforms;
124
125
126
127
128
    if (platform == NULL) {
        char* defaultPlatform = getenv("OPENMM_DEFAULT_PLATFORM");
        if (defaultPlatform != NULL)
            platform = &Platform::getPlatformByName(string(defaultPlatform));
    }
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
    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;
147
            platform->contextCreated(*this, validatedProperties);
148
149
150
151
152
153
154
155
            break;
        }
        catch (...) {
            if (i > 0)
                continue;
            throw;
        }
    }
156
157
158
    
    // Create and initialize kernels and other objects.
    
159
    initializeForcesKernel = platform->createKernel(CalcForcesAndEnergyKernel::Name(), *this);
160
    initializeForcesKernel.getAs<CalcForcesAndEnergyKernel>().initialize(system);
161
    updateStateDataKernel = platform->createKernel(UpdateStateDataKernel::Name(), *this);
162
    updateStateDataKernel.getAs<UpdateStateDataKernel>().initialize(system);
163
    applyConstraintsKernel = platform->createKernel(ApplyConstraintsKernel::Name(), *this);
164
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().initialize(system);
165
    virtualSitesKernel = platform->createKernel(VirtualSitesKernel::Name(), *this);
166
    virtualSitesKernel.getAs<VirtualSitesKernel>().initialize(system);
167
168
    Vec3 periodicBoxVectors[3];
    system.getDefaultPeriodicBoxVectors(periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
169
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
170
    for (size_t i = 0; i < forceImpls.size(); ++i)
171
        forceImpls[i]->initialize(*this);
172
    integrator.initialize(*this);
173
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, vector<Vec3>(system.getNumParticles()));
174
175
}

176
ContextImpl::~ContextImpl() {
177
178
    for (int i = 0; i < (int) forceImpls.size(); ++i)
        delete forceImpls[i];
179
180
181
182
183
184
185
    
    // Make sure all kernels get properly deleted before contextDestroyed() is called.
    
    initializeForcesKernel = Kernel();
    updateStateDataKernel = Kernel();
    applyConstraintsKernel = Kernel();
    virtualSitesKernel = Kernel();
186
187
188
189
190
191
    if (!integratorIsDeleted) {
        // The Context is being deleted before the Integrator, so call cleanup() on it now.
        
        integrator.cleanup();
        integrator.context = NULL;
    }
192
    platform->contextDestroyed(*this);
193
194
}

195
double ContextImpl::getTime() const {
196
    return updateStateDataKernel.getAs<const UpdateStateDataKernel>().getTime(*this);
197
198
}

199
void ContextImpl::setTime(double t) {
200
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setTime(*this, t);
201
202
203
}

void ContextImpl::getPositions(std::vector<Vec3>& positions) {
204
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPositions(*this, positions);
205
206
207
}

void ContextImpl::setPositions(const std::vector<Vec3>& positions) {
208
    hasSetPositions = true;
209
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPositions(*this, positions);
210
    integrator.stateChanged(State::Positions);
211
212
213
}

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

void ContextImpl::setVelocities(const std::vector<Vec3>& velocities) {
218
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, velocities);
219
    integrator.stateChanged(State::Velocities);
220
221
222
}

void ContextImpl::getForces(std::vector<Vec3>& forces) {
223
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getForces(*this, forces);
224
225
}

226
227
228
229
const std::map<std::string, double>& ContextImpl::getParameters() const {
    return parameters;
}

230
double ContextImpl::getParameter(std::string name) {
231
    if (parameters.find(name) == parameters.end())
232
        throw OpenMMException("Called getParameter() with invalid parameter name: "+name);
233
234
235
    return parameters[name];
}

236
void ContextImpl::setParameter(std::string name, double value) {
237
    if (parameters.find(name) == parameters.end())
238
        throw OpenMMException("Called setParameter() with invalid parameter name: "+name);
239
    parameters[name] = value;
240
    integrator.stateChanged(State::Parameters);
241
242
}

243
void ContextImpl::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) {
244
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPeriodicBoxVectors(*this, a, b, c);
245
246
247
248
249
}

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.");
250
251
252
253
    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.");
254
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, a, b, c);
255
256
}

257
void ContextImpl::applyConstraints(double tol) {
258
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().apply(*this, tol);
259
260
}

261
262
263
264
void ContextImpl::applyVelocityConstraints(double tol) {
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().applyToVelocities(*this, tol);
}

265
void ContextImpl::computeVirtualSites() {
266
    virtualSitesKernel.getAs<VirtualSitesKernel>().computePositions(*this);
267
268
}

269
double ContextImpl::calcForcesAndEnergy(bool includeForces, bool includeEnergy, int groups) {
270
271
    if (!hasSetPositions)
        throw OpenMMException("Particle positions have not been set");
272
    lastForceGroups = groups;
273
    CalcForcesAndEnergyKernel& kernel = initializeForcesKernel.getAs<CalcForcesAndEnergyKernel>();
274
275
276
277
278
279
280
281
282
283
    while (true) {
        double energy = 0.0;
        kernel.beginComputation(*this, includeForces, includeEnergy, groups);
        for (int i = 0; i < (int) forceImpls.size(); ++i)
            energy += forceImpls[i]->calcForcesAndEnergy(*this, includeForces, includeEnergy, groups);
        bool valid = true;
        energy += kernel.finishComputation(*this, includeForces, includeEnergy, groups, valid);
        if (valid)
            return energy;
    }
284
285
}

286
287
288
289
int ContextImpl::getLastForceGroups() const {
    return lastForceGroups;
}

290
double ContextImpl::calcKineticEnergy() {
291
    return integrator.computeKineticEnergy();
292
293
}

294
void ContextImpl::updateContextState() {
295
296
297
298
    for (int i = 0; i < (int) forceImpls.size(); ++i)
        forceImpls[i]->updateContextState(*this);
}

299
300
301
302
const vector<ForceImpl*>& ContextImpl::getForceImpls() const {
    return forceImpls;
}

303
304
305
306
vector<ForceImpl*>& ContextImpl::getForceImpls() {
    return forceImpls;
}

307
void* ContextImpl::getPlatformData() {
308
309
310
    return platformData;
}

311
312
313
314
const void* ContextImpl::getPlatformData() const {
    return platformData;
}

315
void ContextImpl::setPlatformData(void* data) {
316
317
    platformData = data;
}
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337

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));
    }
    for (int i = 0; i < (int) forceImpls.size(); i++) {
        vector<pair<int, int> > forceBonds = forceImpls[i]->getBondedParticles();
        bonds.insert(bonds.end(), forceBonds.begin(), forceBonds.end());
    }
338
339
340
341
342
343
344
    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)));
        }
    }
345
346
347
348
349
350
351
352
353
354

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

    int numParticles = system.getNumParticles();
    vector<vector<int> > particleBonds(numParticles);
    for (int i = 0; i < (int) bonds.size(); i++) {
        particleBonds[bonds[i].first].push_back(bonds[i].second);
        particleBonds[bonds[i].second].push_back(bonds[i].first);
    }

355
    // Now identify particles by which molecule they belong to.
356

357
358
359
360
361
362
363
364
365
    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.
    
366
367
368
    vector<int> particleMolecule(numParticles, -1);
    int numMolecules = 0;
    for (int i = 0; i < numParticles; i++)
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
        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);
400
401
402
403
404
    for (int i = 0; i < numParticles; i++)
        molecules[particleMolecule[i]].push_back(i);
    return molecules;
}

Peter Eastman's avatar
Peter Eastman committed
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
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
420
    stream.write(CHECKPOINT_MAGIC_BYTES, sizeof(CHECKPOINT_MAGIC_BYTES)/sizeof(CHECKPOINT_MAGIC_BYTES[0]));
Peter Eastman's avatar
Peter Eastman committed
421
422
423
424
425
426
427
428
429
    writeString(stream, getPlatform().getName());
    int numParticles = getSystem().getNumParticles();
    stream.write((char*) &numParticles, sizeof(int));
    int numParameters = parameters.size();
    stream.write((char*) &numParameters, sizeof(int));
    for (map<string, double>::const_iterator iter = parameters.begin(); iter != parameters.end(); ++iter) {
        writeString(stream, iter->first);
        stream.write((char*) &iter->second, sizeof(double));
    }
430
    updateStateDataKernel.getAs<UpdateStateDataKernel>().createCheckpoint(*this, stream);
Peter Eastman's avatar
Peter Eastman committed
431
432
433
434
    stream.flush();
}

void ContextImpl::loadCheckpoint(istream& stream) {
Robert McGibbon's avatar
Robert McGibbon committed
435
436
437
438
    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
439
440
        throw OpenMMException("loadCheckpoint: Checkpoint header was not correct");

Peter Eastman's avatar
Peter Eastman committed
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
    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;
    }
456
    updateStateDataKernel.getAs<UpdateStateDataKernel>().loadCheckpoint(*this, stream);
457
    hasSetPositions = true;
Peter Eastman's avatar
Peter Eastman committed
458
}