ContextImpl.cpp 19.2 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-2013 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
87
88
89
90
91
92
93
94
95
96
97
    // Validate the list of properties.

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

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

190
double ContextImpl::getTime() const {
191
    return updateStateDataKernel.getAs<const UpdateStateDataKernel>().getTime(*this);
192
193
}

194
void ContextImpl::setTime(double t) {
195
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setTime(*this, t);
196
197
198
}

void ContextImpl::getPositions(std::vector<Vec3>& positions) {
199
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPositions(*this, positions);
200
201
202
}

void ContextImpl::setPositions(const std::vector<Vec3>& positions) {
203
    hasSetPositions = true;
204
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPositions(*this, positions);
205
    integrator.stateChanged(State::Positions);
206
207
208
}

void ContextImpl::getVelocities(std::vector<Vec3>& velocities) {
209
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getVelocities(*this, velocities);
210
211
212
}

void ContextImpl::setVelocities(const std::vector<Vec3>& velocities) {
213
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, velocities);
214
    integrator.stateChanged(State::Velocities);
215
216
217
}

void ContextImpl::getForces(std::vector<Vec3>& forces) {
218
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getForces(*this, forces);
219
220
}

221
222
223
224
const std::map<std::string, double>& ContextImpl::getParameters() const {
    return parameters;
}

225
double ContextImpl::getParameter(std::string name) {
226
    if (parameters.find(name) == parameters.end())
227
        throw OpenMMException("Called getParameter() with invalid parameter name: "+name);
228
229
230
    return parameters[name];
}

231
void ContextImpl::setParameter(std::string name, double value) {
232
    if (parameters.find(name) == parameters.end())
233
        throw OpenMMException("Called setParameter() with invalid parameter name: "+name);
234
    parameters[name] = value;
235
    integrator.stateChanged(State::Parameters);
236
237
}

238
void ContextImpl::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) {
239
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPeriodicBoxVectors(*this, a, b, c);
240
241
242
243
244
}

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.");
245
246
247
248
    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.");
249
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, a, b, c);
250
251
}

252
void ContextImpl::applyConstraints(double tol) {
253
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().apply(*this, tol);
254
255
}

256
257
258
259
void ContextImpl::applyVelocityConstraints(double tol) {
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().applyToVelocities(*this, tol);
}

260
void ContextImpl::computeVirtualSites() {
261
    virtualSitesKernel.getAs<VirtualSitesKernel>().computePositions(*this);
262
263
}

264
double ContextImpl::calcForcesAndEnergy(bool includeForces, bool includeEnergy, int groups) {
265
266
    if (!hasSetPositions)
        throw OpenMMException("Particle positions have not been set");
267
    lastForceGroups = groups;
268
    CalcForcesAndEnergyKernel& kernel = initializeForcesKernel.getAs<CalcForcesAndEnergyKernel>();
269
270
271
272
273
274
275
276
277
278
    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;
    }
279
280
}

281
282
283
284
int ContextImpl::getLastForceGroups() const {
    return lastForceGroups;
}

285
double ContextImpl::calcKineticEnergy() {
286
    return integrator.computeKineticEnergy();
287
288
}

289
void ContextImpl::updateContextState() {
290
291
292
293
    for (int i = 0; i < (int) forceImpls.size(); ++i)
        forceImpls[i]->updateContextState(*this);
}

294
295
296
297
const vector<ForceImpl*>& ContextImpl::getForceImpls() const {
    return forceImpls;
}

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

302
void* ContextImpl::getPlatformData() {
303
304
305
    return platformData;
}

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

310
void ContextImpl::setPlatformData(void* data) {
311
312
    platformData = data;
}
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332

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());
    }
333
334
335
336
337
338
339
    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)));
        }
    }
340
341
342
343
344
345
346
347
348
349

    // 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);
    }

350
    // Now identify particles by which molecule they belong to.
351

352
353
354
355
356
357
358
359
360
    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.
    
361
362
363
    vector<int> particleMolecule(numParticles, -1);
    int numMolecules = 0;
    for (int i = 0; i < numParticles; i++)
364
365
366
367
368
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
        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);
395
396
397
398
399
    for (int i = 0; i < numParticles; i++)
        molecules[particleMolecule[i]].push_back(i);
    return molecules;
}

Peter Eastman's avatar
Peter Eastman committed
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
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
415
    stream.write(CHECKPOINT_MAGIC_BYTES, sizeof(CHECKPOINT_MAGIC_BYTES)/sizeof(CHECKPOINT_MAGIC_BYTES[0]));
Peter Eastman's avatar
Peter Eastman committed
416
417
418
419
420
421
422
423
424
    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));
    }
425
    updateStateDataKernel.getAs<UpdateStateDataKernel>().createCheckpoint(*this, stream);
Peter Eastman's avatar
Peter Eastman committed
426
427
428
429
    stream.flush();
}

void ContextImpl::loadCheckpoint(istream& stream) {
Robert McGibbon's avatar
Robert McGibbon committed
430
431
432
433
    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
434
435
        throw OpenMMException("loadCheckpoint: Checkpoint header was not correct");

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