ReferenceCustomDynamics.cpp 16.4 KB
Newer Older
1

2
/* Portions copyright (c) 2011-2013 Stanford University and Simbios.
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
 * Contributors: Peter Eastman
 *
 * 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.
 */

25
26
#include "SimTKOpenMMCommon.h"
#include "SimTKOpenMMUtilities.h"
27
#include "ReferenceVirtualSites.h"
28
#include "ReferenceCustomDynamics.h"
29
#include "openmm/OpenMMException.h"
30
31
32
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/ForceImpl.h"
#include "lepton/Operation.h"
33
34
#include "lepton/ParsedExpression.h"
#include "lepton/Parser.h"
35
#include <set>
36
#include <sstream>
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

using namespace std;
using namespace OpenMM;

/**---------------------------------------------------------------------------------------

   ReferenceCustomDynamics constructor

   @param numberOfAtoms  number of atoms
   @param integrator     the integrator definition to use

   --------------------------------------------------------------------------------------- */

ReferenceCustomDynamics::ReferenceCustomDynamics(int numberOfAtoms, const CustomIntegrator& integrator) : 
           ReferenceDynamics(numberOfAtoms, integrator.getStepSize(), 0.0), integrator(integrator) {
52
    sumBuffer.resize(numberOfAtoms);
53
    oldPos.resize(numberOfAtoms);
54
55
56
57
58
59
60
    stepType.resize(integrator.getNumComputations());
    stepVariable.resize(integrator.getNumComputations());
    stepExpression.resize(integrator.getNumComputations());
    for (int i = 0; i < integrator.getNumComputations(); i++) {
        string expression;
        integrator.getComputationStep(i, stepType[i], stepVariable[i], expression);
        if (expression.length() > 0)
61
62
63
64
65
66
67
68
            stepExpression[i] = Lepton::Parser::parse(expression).optimize().createProgram();
    }
    kineticEnergyExpression = Lepton::Parser::parse(integrator.getKineticEnergyExpression()).optimize().createProgram();
    kineticEnergyNeedsForce = false;
    for (int i = 0; i < kineticEnergyExpression.getNumOperations(); i++) {
        const Lepton::Operation& op = kineticEnergyExpression.getOperation(i);
        if (op.getId() == Lepton::Operation::VARIABLE && op.getName() == "f")
            kineticEnergyNeedsForce = true;
69
    }
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
}

/**---------------------------------------------------------------------------------------

   ReferenceCustomDynamics destructor

   --------------------------------------------------------------------------------------- */

ReferenceCustomDynamics::~ReferenceCustomDynamics() {
}

/**---------------------------------------------------------------------------------------

   Update -- driver routine for performing Custom dynamics update of coordinates
   and velocities

   @param context             the context this integrator is updating
   @param numberOfAtoms       number of atoms
   @param atomCoordinates     atom coordinates
   @param velocities          velocities
   @param forces              forces
   @param masses              atom masses
   @param globals             a map containing values of global variables
   @param forcesAreValid      whether the current forces are valid or need to be recomputed

   --------------------------------------------------------------------------------------- */

void ReferenceCustomDynamics::update(ContextImpl& context, int numberOfAtoms, vector<RealVec>& atomCoordinates,
                                     vector<RealVec>& velocities, vector<RealVec>& forces, vector<RealOpenMM>& masses,
99
                                     map<string, RealOpenMM>& globals, vector<vector<RealVec> >& perDof, bool& forcesAreValid, RealOpenMM tolerance){
100
101
    int numSteps = stepType.size();
    globals.insert(context.getParameters().begin(), context.getParameters().end());
102
    oldPos = atomCoordinates;
103
104
105
106
107
108
109
    if (invalidatesForces.size() == 0) {
        // The first time this is called, work out when to recompute forces and energy.  First build a
        // list of every step that invalidates the forces.
        
        invalidatesForces.resize(numSteps, false);
        needsForces.resize(numSteps, false);
        needsEnergy.resize(numSteps, false);
110
111
        forceGroup.resize(numSteps, -2);
        forceName.resize(numSteps, "f");
112
        energyName.resize(numSteps, "energy");
113
114
115
116
117
118
119
120
121
122
123
124
        set<string> affectsForce;
        affectsForce.insert("x");
        for (vector<ForceImpl*>::const_iterator iter = context.getForceImpls().begin(); iter != context.getForceImpls().end(); ++iter) {
            const map<string, double> params = (*iter)->getDefaultParameters();
            for (map<string, double>::const_iterator param = params.begin(); param != params.end(); ++param)
                affectsForce.insert(param->first);
        }
        for (int i = 0; i < numSteps; i++)
            invalidatesForces[i] = (stepType[i] == CustomIntegrator::ConstrainPositions || affectsForce.find(stepVariable[i]) != affectsForce.end());
        
        // Make a list of which steps require valid forces or energy to be known.
        
125
        vector<string> forceGroupName;
126
        vector<string> energyGroupName;
127
        for (int i = 0; i < 32; i++) {
128
129
130
131
132
133
            stringstream fname;
            fname << "f" << i;
            forceGroupName.push_back(fname.str());
            stringstream ename;
            ename << "energy" << i;
            energyGroupName.push_back(ename.str());
134
        }
135
136
137
138
139
        for (int i = 0; i < numSteps; i++) {
            if (stepType[i] == CustomIntegrator::ComputeGlobal || stepType[i] == CustomIntegrator::ComputePerDof || stepType[i] == CustomIntegrator::ComputeSum) {
                for (int j = 0; j < stepExpression[i].getNumOperations(); j++) {
                    const Lepton::Operation& op = stepExpression[i].getOperation(j);
                    if (op.getId() == Lepton::Operation::VARIABLE) {
140
141
142
                        if (op.getName() == "energy") {
                            if (forceGroup[i] != -2)
                                throw OpenMMException("A single computation step cannot depend on multiple force groups");
143
                            needsEnergy[i] = true;
144
145
                            forceGroup[i] = -1;
                        }
146
147
148
149
150
151
152
153
154
155
156
                        else if (op.getName().substr(0, 6) == "energy") {
                            for (int k = 0; k < (int) energyGroupName.size(); k++)
                                if (op.getName() == energyGroupName[k]) {
                                    if (forceGroup[i] != -2)
                                        throw OpenMMException("A single computation step cannot depend on multiple force groups");
                                    needsForces[i] = true;
                                    forceGroup[i] = 1<<k;
                                    energyName[i] = energyGroupName[k];
                                    break;
                                }
                        }
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
                        else if (op.getName() == "f") {
                            if (forceGroup[i] != -2)
                                throw OpenMMException("A single computation step cannot depend on multiple force groups");
                            needsForces[i] = true;
                            forceGroup[i] = -1;
                        }
                        else if (op.getName()[0] == 'f') {
                            for (int k = 0; k < (int) forceGroupName.size(); k++)
                                if (op.getName() == forceGroupName[k]) {
                                    if (forceGroup[i] != -2)
                                        throw OpenMMException("A single computation step cannot depend on multiple force groups");
                                    needsForces[i] = true;
                                    forceGroup[i] = 1<<k;
                                    forceName[i] = forceGroupName[k];
                                    break;
                                }
                        }
174
175
176
177
178
179
180
181
                    }
                }
            }
        }
        
        // Build the list of inverse masses.
        
        inverseMasses.resize(numberOfAtoms);
182
183
184
185
186
187
        for (int i = 0; i < numberOfAtoms; i++) {
            if (masses[i] == 0.0)
                inverseMasses[i] = 0.0;
            else
                inverseMasses[i] = 1.0/masses[i];
        }
188
189
190
191
192
    }
    
    // Loop over steps and execute them.
    
    for (int i = 0; i < numSteps; i++) {
193
        if ((needsForces[i] || needsEnergy[i]) && (!forcesAreValid || context.getLastForceGroups() != forceGroup[i])) {
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
            // Recompute forces and or energy.  Figure out what is actually needed
            // between now and the next time they get invalidated again.
            
            bool computeForce = false, computeEnergy = false;
            for (int j = i; ; j++) {
                if (needsForces[j])
                    computeForce = true;
                if (needsEnergy[j])
                    computeEnergy = true;
                if (invalidatesForces[j])
                    break;
                if (j == numSteps-1)
                    j = -1;
                if (j == i-1)
                    break;
            }
            recordChangedParameters(context, globals);
211
            RealOpenMM e = context.calcForcesAndEnergy(computeForce, computeEnergy, forceGroup[i]);
212
213
214
215
            if (computeEnergy)
                energy = e;
            forcesAreValid = true;
        }
216
        globals[energyName[i]] = energy;
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
        
        // Execute the step.
        
        switch (stepType[i]) {
            case CustomIntegrator::ComputeGlobal: {
                map<string, RealOpenMM> variables = globals;
                variables["uniform"] = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber();
                variables["gaussian"] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
                globals[stepVariable[i]] = stepExpression[i].evaluate(variables);
                break;
            }
            case CustomIntegrator::ComputePerDof: {
                vector<RealVec>* results = NULL;
                if (stepVariable[i] == "x")
                    results = &atomCoordinates;
                else if (stepVariable[i] == "v")
                    results = &velocities;
                else {
                    for (int j = 0; j < integrator.getNumPerDofVariables(); j++)
                        if (stepVariable[i] == integrator.getPerDofVariableName(j))
                            results = &perDof[j];
                }
                if (results == NULL)
                    throw OpenMMException("Illegal per-DOF output variable: "+stepVariable[i]);
241
                computePerDof(numberOfAtoms, *results, atomCoordinates, velocities, forces, masses, globals, perDof, stepExpression[i], forceName[i]);
242
243
244
                break;
            }
            case CustomIntegrator::ComputeSum: {
245
                computePerDof(numberOfAtoms, sumBuffer, atomCoordinates, velocities, forces, masses, globals, perDof, stepExpression[i], forceName[i]);
246
247
                RealOpenMM sum = 0.0;
                for (int j = 0; j < numberOfAtoms; j++)
248
249
                    if (masses[j] != 0.0)
                        sum += sumBuffer[j][0]+sumBuffer[j][1]+sumBuffer[j][2];
250
251
252
253
                globals[stepVariable[i]] = sum;
                break;
            }
            case CustomIntegrator::ConstrainPositions: {
254
                getReferenceConstraintAlgorithm()->apply(oldPos, atomCoordinates, inverseMasses, tolerance);
255
                oldPos = atomCoordinates;
256
257
258
                break;
            }
            case CustomIntegrator::ConstrainVelocities: {
259
                getReferenceConstraintAlgorithm()->applyToVelocities(oldPos, velocities, inverseMasses, tolerance);
260
                break;
261
262
263
264
265
266
267
268
269
270
            }
            case CustomIntegrator::UpdateContextState: {
                recordChangedParameters(context, globals);
                context.updateContextState();
                globals.insert(context.getParameters().begin(), context.getParameters().end());
            }
        }
        if (invalidatesForces[i])
            forcesAreValid = false;
    }
271
    ReferenceVirtualSites::computePositions(context.getSystem(), atomCoordinates);
272
273
274
275
    incrementTimeStep();
    recordChangedParameters(context, globals);
}

276
void ReferenceCustomDynamics::computePerDof(int numberOfAtoms, vector<RealVec>& results, const vector<RealVec>& atomCoordinates,
277
              const vector<RealVec>& velocities, const vector<RealVec>& forces, const vector<RealOpenMM>& masses,
278
279
              const map<string, RealOpenMM>& globals, const vector<vector<RealVec> >& perDof,
              const Lepton::ExpressionProgram& expression, const std::string& forceName) {
280
281
282
283
    // Loop over all degrees of freedom.
    
    map<string, RealOpenMM> variables = globals;
    for (int i = 0; i < numberOfAtoms; i++) {
284
285
286
287
288
289
290
        if (masses[i] != 0.0) {
            variables["m"] = masses[i];
            for (int j = 0; j < 3; j++) {
                // Compute the expression.

                variables["x"] = atomCoordinates[i][j];
                variables["v"] = velocities[i][j];
291
                variables[forceName] = forces[i][j];
292
293
294
295
296
297
                variables["uniform"] = SimTKOpenMMUtilities::getUniformlyDistributedRandomNumber();
                variables["gaussian"] = SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
                for (int k = 0; k < (int) perDof.size(); k++)
                    variables[integrator.getPerDofVariableName(k)] = perDof[k][i][j];
                results[i][j] = expression.evaluate(variables);
            }
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
        }
    }
}

/**
 * Check which context parameters have changed and register them with the context.
 */
void ReferenceCustomDynamics::recordChangedParameters(OpenMM::ContextImpl& context, std::map<std::string, RealOpenMM>& globals) {
    for (map<string, double>::const_iterator iter = context.getParameters().begin(); iter != context.getParameters().end(); ++iter) {
        string name = iter->first;
        double value = globals[name];
        if (value != iter->second)
            context.setParameter(name, globals[name]);
    }
}
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344

/**---------------------------------------------------------------------------------------

   Compute the kinetic energy of the system.

   @param context             the context this integrator is updating
   @param numberOfAtoms       number of atoms
   @param atomCoordinates     atom coordinates
   @param velocities          velocities
   @param forces              forces
   @param masses              atom masses
   @param globals             a map containing values of global variables
   @param perDof              the values of per-DOF variables
   @param forcesAreValid      whether the current forces are valid or need to be recomputed

   --------------------------------------------------------------------------------------- */

double ReferenceCustomDynamics::computeKineticEnergy(OpenMM::ContextImpl& context, int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
        std::vector<OpenMM::RealVec>& velocities, std::vector<OpenMM::RealVec>& forces, std::vector<RealOpenMM>& masses,
        std::map<std::string, RealOpenMM>& globals, std::vector<std::vector<OpenMM::RealVec> >& perDof, bool& forcesAreValid) {
    globals.insert(context.getParameters().begin(), context.getParameters().end());
    if (kineticEnergyNeedsForce) {
        energy = context.calcForcesAndEnergy(true, true, -1);
        forcesAreValid = true;
    }
    computePerDof(numberOfAtoms, sumBuffer, atomCoordinates, velocities, forces, masses, globals, perDof, kineticEnergyExpression, "f");
    RealOpenMM sum = 0.0;
    for (int j = 0; j < numberOfAtoms; j++)
        if (masses[j] != 0.0)
            sum += sumBuffer[j][0]+sumBuffer[j][1]+sumBuffer[j][2];
    return sum;
}