Context.cpp 9.88 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
#include "openmm/Context.h"
#include "openmm/internal/ContextImpl.h"
34
#include "openmm/OpenMMException.h"
35
#include "openmm/internal/ForceImpl.h"
36
#include "SimTKOpenMMRealType.h"
37
#include "sfmt/SFMT.h"
38
#include <cmath>
39
40
41
42

using namespace OpenMM;
using namespace std;

43
Context::Context(const System& system, Integrator& integrator) : properties(map<string, string>()) {
44
    impl = new ContextImpl(*this, system, integrator, 0, properties);
45
    impl->initialize();
46
47
}

48
Context::Context(const System& system, Integrator& integrator, Platform& platform) : properties(map<string, string>()) {
49
    impl = new ContextImpl(*this, system, integrator, &platform, properties);
50
    impl->initialize();
51
52
}

53
Context::Context(const System& system, Integrator& integrator, Platform& platform, const map<string, string>& properties) : properties(properties) {
54
    impl = new ContextImpl(*this, system, integrator, &platform, properties);
55
    impl->initialize();
56
57
}

58
Context::~Context() {
Lee-Ping's avatar
Lee-Ping committed
59
    delete impl;
60
61
}

62
const System& Context::getSystem() const {
63
64
65
66
    return impl->getSystem();

}

67
const Integrator& Context::getIntegrator() const {
68
69
70
    return impl->getIntegrator();
}

71
Integrator& Context::getIntegrator() {
72
    return impl->getIntegrator();
73
74
}

75
const Platform& Context::getPlatform() const {
76
77
    return impl->getPlatform();
}
78

79
Platform& Context::getPlatform() {
80
    return impl->getPlatform();
81
82
}

83
State Context::getState(int types, bool enforcePeriodicBox, int groups) const {
84
    State::StateBuilder builder(impl->getTime());
85
86
    Vec3 periodicBoxSize[3];
    impl->getPeriodicBoxVectors(periodicBoxSize[0], periodicBoxSize[1], periodicBoxSize[2]);
87
    builder.setPeriodicBoxVectors(periodicBoxSize[0], periodicBoxSize[1], periodicBoxSize[2]);
88
89
    bool includeForces = types&State::Forces;
    bool includeEnergy = types&State::Energy;
90
91
92
    bool includeParameterDerivs = types&State::ParameterDerivatives;
    if (includeForces || includeEnergy || includeParameterDerivs) {
        double energy = impl->calcForcesAndEnergy(includeForces || includeEnergy || includeParameterDerivs, includeEnergy, groups);
93
        if (includeEnergy)
94
95
96
97
98
99
            builder.setEnergy(impl->calcKineticEnergy(), energy);
        if (includeForces) {
            vector<Vec3> forces;
            impl->getForces(forces);
            builder.setForces(forces);
        }
100
    }
101
    if (types&State::Parameters) {
102
        map<string, double> params;
peastman's avatar
peastman committed
103
104
        for (auto& param : impl->parameters)
            params[param.first] = param.second;
105
        builder.setParameters(params);
106
    }
107
108
109
110
111
    if (types&State::ParameterDerivatives) {
        map<string, double> derivs;
        impl->getEnergyParameterDerivatives(derivs);
        builder.setEnergyParameterDerivatives(derivs);
    }
112
    if (types&State::Positions) {
113
114
        vector<Vec3> positions;
        impl->getPositions(positions);
115
116
        if (enforcePeriodicBox) {
            const vector<vector<int> >& molecules = impl->getMolecules();
peastman's avatar
peastman committed
117
            for (auto& mol : molecules) {
118
119
120
                // Find the molecule center.

                Vec3 center;
peastman's avatar
peastman committed
121
122
123
                for (int j : mol)
                    center += positions[j];
                center *= 1.0/mol.size();
124
125

                // Find the displacement to move it into the first periodic box.
126
                Vec3 diff;
127
128
129
                diff += periodicBoxSize[2]*floor(center[2]/periodicBoxSize[2][2]);
                diff += periodicBoxSize[1]*floor((center[1]-diff[1])/periodicBoxSize[1][1]);
                diff += periodicBoxSize[0]*floor((center[0]-diff[0])/periodicBoxSize[0][0]);
130
131

                // Translate all the particles in the molecule.
peastman's avatar
peastman committed
132
133
                for (int j : mol)
                    positions[j] -= diff;
134
135
            }
        }
136
137
138
139
140
141
        builder.setPositions(positions);
    }
    if (types&State::Velocities) {
        vector<Vec3> velocities;
        impl->getVelocities(velocities);
        builder.setVelocities(velocities);
142
    }
143
    return builder.getState();
144
145
}

Peter Eastman's avatar
Peter Eastman committed
146
147
148
149
150
void Context::setState(const State& state) {
    setTime(state.getTime());
    Vec3 a, b, c;
    state.getPeriodicBoxVectors(a, b, c);
    setPeriodicBoxVectors(a, b, c);
peastman's avatar
peastman committed
151
    if ((state.getDataTypes()&State::Positions) != 0)
Peter Eastman's avatar
Peter Eastman committed
152
        setPositions(state.getPositions());
peastman's avatar
peastman committed
153
    if ((state.getDataTypes()&State::Velocities) != 0)
Peter Eastman's avatar
Peter Eastman committed
154
        setVelocities(state.getVelocities());
peastman's avatar
peastman committed
155
    if ((state.getDataTypes()&State::Parameters) != 0)
peastman's avatar
peastman committed
156
157
        for (auto& param : state.getParameters())
            setParameter(param.first, param.second);
Peter Eastman's avatar
Peter Eastman committed
158
159
}

160
void Context::setTime(double time) {
161
162
163
    impl->setTime(time);
}

164
void Context::setPositions(const vector<Vec3>& positions) {
Peter Eastman's avatar
Peter Eastman committed
165
    if ((int) positions.size() != impl->getSystem().getNumParticles())
166
        throw OpenMMException("Called setPositions() on a Context with the wrong number of positions");
167
    impl->setPositions(positions);
168
169
}

170
void Context::setVelocities(const vector<Vec3>& velocities) {
Peter Eastman's avatar
Peter Eastman committed
171
    if ((int) velocities.size() != impl->getSystem().getNumParticles())
172
        throw OpenMMException("Called setVelocities() on a Context with the wrong number of velocities");
173
    impl->setVelocities(velocities);
174
175
}

176
void Context::setVelocitiesToTemperature(double temperature, int randomSeed) {
177
    const System& system = impl->getSystem();
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
    
    // Generate the list of Gaussian random numbers.
    
    OpenMM_SFMT::SFMT sfmt;
    init_gen_rand(randomSeed, sfmt);
    vector<double> randoms;
    while (randoms.size() < system.getNumParticles()*3) {
        double x, y, r2;
        do {
            x = 2.0*genrand_real2(sfmt)-1.0;
            y = 2.0*genrand_real2(sfmt)-1.0;
            r2 = x*x + y*y;
        } while (r2 >= 1.0 || r2 == 0.0);
        double multiplier = sqrt((-2.0*log(r2))/r2);
        randoms.push_back(x*multiplier);
        randoms.push_back(y*multiplier);
    }
    
    // Assign the velocities.
    
    vector<Vec3> velocities(system.getNumParticles(), Vec3());
    int nextRandom = 0;
    for (int i = 0; i < system.getNumParticles(); i++) {
        double mass = system.getParticleMass(i);
        if (mass != 0) {
            double velocityScale = sqrt(BOLTZ*temperature/mass);
            velocities[i] = Vec3(randoms[nextRandom++], randoms[nextRandom++], randoms[nextRandom++])*velocityScale;
        }
    }
    setVelocities(velocities);
    impl->applyVelocityConstraints(1e-5);
}

Peter Eastman's avatar
Peter Eastman committed
211
212
213
214
const map<string, double>& Context::getParameters() const {
    return impl->getParameters();
}

215
double Context::getParameter(const string& name) const {
216
217
218
    return impl->getParameter(name);
}

219
void Context::setParameter(const string& name, double value) {
220
221
222
    impl->setParameter(name, value);
}

223
224
void Context::setPeriodicBoxVectors(const Vec3& a, const Vec3& b, const Vec3& c) {
    impl->setPeriodicBoxVectors(a, b, c);
225
226
}

227
228
229
230
void Context::applyConstraints(double tol) {
    impl->applyConstraints(tol);
}

231
232
233
234
void Context::applyVelocityConstraints(double tol) {
    impl->applyVelocityConstraints(tol);
}

235
236
237
238
void Context::computeVirtualSites() {
    impl->computeVirtualSites();
}

239
void Context::reinitialize() {
240
    const System& system = impl->getSystem();
241
242
    Integrator& integrator = impl->getIntegrator();
    Platform& platform = impl->getPlatform();
243
    integrator.cleanup();
244
    delete impl;
245
    impl = new ContextImpl(*this, system, integrator, &platform, properties);
246
    impl->initialize();
247
}
Peter Eastman's avatar
Peter Eastman committed
248
249
250
251
252
253
254
255

void Context::createCheckpoint(ostream& stream) {
    impl->createCheckpoint(stream);
}

void Context::loadCheckpoint(istream& stream) {
    impl->loadCheckpoint(stream);
}
256
257
258
259

ContextImpl& Context::getImpl() {
    return *impl;
}
260

261
262
263
264
const ContextImpl& Context::getImpl() const {
    return *impl;
}

265
266
267
const vector<vector<int> >& Context::getMolecules() const {
    return impl->getMolecules();
}