ContextImpl.cpp 15.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-2012 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
42
#include "openmm/Context.h"
#include <iostream>
43
#include <map>
44
#include <utility>
45
46
47
#include <vector>

using namespace OpenMM;
Peter Eastman's avatar
Peter Eastman committed
48
using namespace std;
49

50
ContextImpl::ContextImpl(Context& owner, System& system, Integrator& integrator, Platform* platform, const map<string, string>& properties) :
51
         owner(owner), system(system), integrator(integrator), hasInitializedForces(false), lastForceGroups(-1), platform(platform), platformData(NULL) {
52
53
    if (system.getNumParticles() == 0)
        throw OpenMMException("Cannot create a Context for a System with no particles");
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
    
    // 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);
        if (system.getParticleMass(particle1) == 0.0 || system.getParticleMass(particle2) == 0.0)
            throw OpenMMException("A constraint cannot involve a massless particle");
    }
    
    // Find the list of kernels required.
    
77
78
    vector<string> kernelNames;
    kernelNames.push_back(CalcKineticEnergyKernel::Name());
79
    kernelNames.push_back(CalcForcesAndEnergyKernel::Name());
80
    kernelNames.push_back(UpdateStateDataKernel::Name());
81
    for (int i = 0; i < system.getNumForces(); ++i) {
82
        forceImpls.push_back(system.getForce(i).createImpl());
83
84
        map<string, double> forceParameters = forceImpls[forceImpls.size()-1]->getDefaultParameters();
        parameters.insert(forceParameters.begin(), forceParameters.end());
85
86
87
        vector<string> forceKernels = forceImpls[forceImpls.size()-1]->getKernelNames();
        kernelNames.insert(kernelNames.begin(), forceKernels.begin(), forceKernels.end());
    }
88
    hasInitializedForces = true;
89
90
    vector<string> integratorKernels = integrator.getKernelNames();
    kernelNames.insert(kernelNames.begin(), integratorKernels.begin(), integratorKernels.end());
91
    if (platform == 0)
92
        this->platform = platform = &Platform::findPlatform(kernelNames);
93
    else if (!platform->supportsKernels(kernelNames))
94
        throw OpenMMException("Specified a Platform for a Context which does not support all required kernels");
95
96
97
    
    // Create and initialize kernels and other objects.
    
98
    platform->contextCreated(*this, properties);
99
    initializeForcesKernel = platform->createKernel(CalcForcesAndEnergyKernel::Name(), *this);
100
    initializeForcesKernel.getAs<CalcForcesAndEnergyKernel>().initialize(system);
101
    kineticEnergyKernel = platform->createKernel(CalcKineticEnergyKernel::Name(), *this);
102
    kineticEnergyKernel.getAs<CalcKineticEnergyKernel>().initialize(system);
103
    updateStateDataKernel = platform->createKernel(UpdateStateDataKernel::Name(), *this);
104
    updateStateDataKernel.getAs<UpdateStateDataKernel>().initialize(system);
105
    applyConstraintsKernel = platform->createKernel(ApplyConstraintsKernel::Name(), *this);
106
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().initialize(system);
107
    virtualSitesKernel = platform->createKernel(VirtualSitesKernel::Name(), *this);
108
    virtualSitesKernel.getAs<VirtualSitesKernel>().initialize(system);
109
110
    Vec3 periodicBoxVectors[3];
    system.getDefaultPeriodicBoxVectors(periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
111
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, periodicBoxVectors[0], periodicBoxVectors[1], periodicBoxVectors[2]);
112
    for (size_t i = 0; i < forceImpls.size(); ++i)
113
        forceImpls[i]->initialize(*this);
114
    integrator.initialize(*this);
115
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, vector<Vec3>(system.getNumParticles()));
116
117
}

118
ContextImpl::~ContextImpl() {
119
120
    for (int i = 0; i < (int) forceImpls.size(); ++i)
        delete forceImpls[i];
121
122
123
124
125
126
127
128
    
    // Make sure all kernels get properly deleted before contextDestroyed() is called.
    
    initializeForcesKernel = Kernel();
    kineticEnergyKernel = Kernel();
    updateStateDataKernel = Kernel();
    applyConstraintsKernel = Kernel();
    virtualSitesKernel = Kernel();
129
    integrator.cleanup();
130
    platform->contextDestroyed(*this);
131
132
}

133
double ContextImpl::getTime() const {
134
    return updateStateDataKernel.getAs<const UpdateStateDataKernel>().getTime(*this);
135
136
}

137
void ContextImpl::setTime(double t) {
138
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setTime(*this, t);
139
140
141
}

void ContextImpl::getPositions(std::vector<Vec3>& positions) {
142
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPositions(*this, positions);
143
144
145
}

void ContextImpl::setPositions(const std::vector<Vec3>& positions) {
146
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPositions(*this, positions);
147
    integrator.stateChanged(State::Positions);
148
149
150
}

void ContextImpl::getVelocities(std::vector<Vec3>& velocities) {
151
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getVelocities(*this, velocities);
152
153
154
}

void ContextImpl::setVelocities(const std::vector<Vec3>& velocities) {
155
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setVelocities(*this, velocities);
156
    integrator.stateChanged(State::Velocities);
157
158
159
}

void ContextImpl::getForces(std::vector<Vec3>& forces) {
160
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getForces(*this, forces);
161
162
}

163
164
165
166
const std::map<std::string, double>& ContextImpl::getParameters() const {
    return parameters;
}

167
double ContextImpl::getParameter(std::string name) {
168
169
170
171
172
    if (parameters.find(name) == parameters.end())
        throw OpenMMException("Called getParameter() with invalid parameter name");
    return parameters[name];
}

173
void ContextImpl::setParameter(std::string name, double value) {
174
175
176
    if (parameters.find(name) == parameters.end())
        throw OpenMMException("Called setParameter() with invalid parameter name");
    parameters[name] = value;
177
    integrator.stateChanged(State::Parameters);
178
179
}

180
void ContextImpl::getPeriodicBoxVectors(Vec3& a, Vec3& b, Vec3& c) {
181
    updateStateDataKernel.getAs<UpdateStateDataKernel>().getPeriodicBoxVectors(*this, a, b, c);
182
183
184
185
186
187
188
189
190
}

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.");
    if (b[0] != 0.0 || b[2] != 0.0)
        throw OpenMMException("Second periodic box vector must be parallel to y.");
    if (c[0] != 0.0 || c[1] != 0.0)
        throw OpenMMException("Third periodic box vector must be parallel to z.");
191
    updateStateDataKernel.getAs<UpdateStateDataKernel>().setPeriodicBoxVectors(*this, a, b, c);
192
193
}

194
void ContextImpl::applyConstraints(double tol) {
195
    applyConstraintsKernel.getAs<ApplyConstraintsKernel>().apply(*this, tol);
196
197
}

198
void ContextImpl::computeVirtualSites() {
199
    virtualSitesKernel.getAs<VirtualSitesKernel>().computePositions(*this);
200
201
}

202
double ContextImpl::calcForcesAndEnergy(bool includeForces, bool includeEnergy, int groups) {
203
    lastForceGroups = groups;
204
    CalcForcesAndEnergyKernel& kernel = initializeForcesKernel.getAs<CalcForcesAndEnergyKernel>();
205
    double energy = 0.0;
206
    kernel.beginComputation(*this, includeForces, includeEnergy, groups);
207
    for (int i = 0; i < (int) forceImpls.size(); ++i)
208
209
        energy += forceImpls[i]->calcForcesAndEnergy(*this, includeForces, includeEnergy, groups);
    energy += kernel.finishComputation(*this, includeForces, includeEnergy, groups);
210
    return energy;
211
212
}

213
214
215
216
int ContextImpl::getLastForceGroups() const {
    return lastForceGroups;
}

217
double ContextImpl::calcKineticEnergy() {
218
    return kineticEnergyKernel.getAs<CalcKineticEnergyKernel>().execute(*this);
219
220
}

221
void ContextImpl::updateContextState() {
222
223
224
225
    for (int i = 0; i < (int) forceImpls.size(); ++i)
        forceImpls[i]->updateContextState(*this);
}

226
227
228
229
const vector<ForceImpl*>& ContextImpl::getForceImpls() const {
    return forceImpls;
}

230
231
232
233
vector<ForceImpl*>& ContextImpl::getForceImpls() {
    return forceImpls;
}

234
void* ContextImpl::getPlatformData() {
235
236
237
    return platformData;
}

238
239
240
241
const void* ContextImpl::getPlatformData() const {
    return platformData;
}

242
void ContextImpl::setPlatformData(void* data) {
243
244
    platformData = data;
}
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264

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());
    }
265
266
267
268
269
270
271
    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)));
        }
    }
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302

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

    // Now tag particles by which molecule they belong to.

    vector<int> particleMolecule(numParticles, -1);
    int numMolecules = 0;
    for (int i = 0; i < numParticles; i++)
        if (particleMolecule[i] == -1)
            tagParticlesInMolecule(i, numMolecules++, particleMolecule, particleBonds);
    molecules.resize(numMolecules);
    for (int i = 0; i < numParticles; i++)
        molecules[particleMolecule[i]].push_back(i);
    return molecules;
}

void ContextImpl::tagParticlesInMolecule(int particle, int molecule, vector<int>& particleMolecule, vector<vector<int> >& particleBonds) {
    // Recursively tag particles as belonging to a particular molecule.

    particleMolecule[particle] = molecule;
    for (int i = 0; i < (int) particleBonds[particle].size(); i++)
        if (particleMolecule[particleBonds[particle][i]] == -1)
            tagParticlesInMolecule(particleBonds[particle][i], molecule, particleMolecule, particleBonds);
}
Peter Eastman's avatar
Peter Eastman committed
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327

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) {
    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));
    }
328
    updateStateDataKernel.getAs<UpdateStateDataKernel>().createCheckpoint(*this, stream);
Peter Eastman's avatar
Peter Eastman committed
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
    stream.flush();
}

void ContextImpl::loadCheckpoint(istream& stream) {
    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;
    }
348
    updateStateDataKernel.getAs<UpdateStateDataKernel>().loadCheckpoint(*this, stream);
Peter Eastman's avatar
Peter Eastman committed
349
}