CpuCustomGBForce.cpp 22.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

/* Portions copyright (c) 2009-2014 Stanford University and Simbios.
 * 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.
 */

#include <string.h>
#include <sstream>

#include "SimTKOpenMMCommon.h"
#include "SimTKOpenMMLog.h"
#include "SimTKOpenMMUtilities.h"
#include "ReferenceForce.h"
#include "CpuCustomGBForce.h"

34
35
using namespace OpenMM;
using namespace std;
36

37
38
39
CpuCustomGBForce::CpuCustomGBForce(int numAtoms, const vector<Lepton::CompiledExpression>& valueExpressions,
                     const vector<vector<Lepton::CompiledExpression> >& valueDerivExpressions,
                     const vector<vector<Lepton::CompiledExpression> >& valueGradientExpressions,
40
                     const vector<string>& valueNames,
41
                     const vector<CustomGBForce::ComputationType>& valueTypes,
42
                     const vector<Lepton::CompiledExpression>& energyExpressions,
43
44
                     const vector<vector<Lepton::CompiledExpression> >& energyDerivExpressions,
                     const vector<vector<Lepton::CompiledExpression> >& energyGradientExpressions,
45
                     const vector<CustomGBForce::ComputationType>& energyTypes,
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
                     const vector<string>& parameterNames) :
            cutoff(false), periodic(false), valueExpressions(valueExpressions), valueDerivExpressions(valueDerivExpressions), valueGradientExpressions(valueGradientExpressions),
            valueNames(valueNames), valueTypes(valueTypes), energyExpressions(energyExpressions), energyDerivExpressions(energyDerivExpressions), energyGradientExpressions(energyGradientExpressions),
            energyTypes(energyTypes), paramNames(parameterNames) {
    for (int i = 0; i < (int) valueExpressions.size(); i++)
        expressionSet.registerExpression(this->valueExpressions[i]);
    for (int i = 0; i < (int) valueDerivExpressions.size(); i++)
        for (int j = 0; j < (int) valueDerivExpressions[i].size(); j++)
            expressionSet.registerExpression(this->valueDerivExpressions[i][j]);
    for (int i = 0; i < (int) valueGradientExpressions.size(); i++)
        for (int j = 0; j < (int) valueGradientExpressions[i].size(); j++)
            expressionSet.registerExpression(this->valueGradientExpressions[i][j]);
    for (int i = 0; i < (int) energyExpressions.size(); i++)
        expressionSet.registerExpression(this->energyExpressions[i]);
    for (int i = 0; i < (int) energyDerivExpressions.size(); i++)
        for (int j = 0; j < (int) energyDerivExpressions[i].size(); j++)
            expressionSet.registerExpression(this->energyDerivExpressions[i][j]);
    for (int i = 0; i < (int) energyGradientExpressions.size(); i++)
        for (int j = 0; j < (int) energyGradientExpressions[i].size(); j++)
            expressionSet.registerExpression(this->energyGradientExpressions[i][j]);
    xindex = expressionSet.getVariableIndex("x");
    yindex = expressionSet.getVariableIndex("y");
    zindex = expressionSet.getVariableIndex("z");
    rindex = expressionSet.getVariableIndex("r");
    for (int i = 0; i < (int) paramNames.size(); i++) {
        paramIndex.push_back(expressionSet.getVariableIndex(paramNames[i]));
        for (int j = 1; j < 3; j++) {
            stringstream name;
            name << paramNames[i] << j;
            particleParamIndex.push_back(expressionSet.getVariableIndex(name.str()));
        }
    }
    for (int i = 0; i < (int) valueNames.size(); i++) {
        valueIndex.push_back(expressionSet.getVariableIndex(valueNames[i]));
        for (int j = 1; j < 3; j++) {
            stringstream name;
            name << valueNames[i] << j;
            particleValueIndex.push_back(expressionSet.getVariableIndex(name.str()));
        }
    }
86
87
88
89
90
91
92
93
94
95
96
    values.resize(valueTypes.size());
    dEdV.resize(valueTypes.size());
    for (int i = 0; i < (int) values.size(); i++) {
        values[i].resize(numAtoms);
        dEdV[i].resize(numAtoms);
    }
    dVdX.resize(valueDerivExpressions.size());
    dVdY.resize(valueDerivExpressions.size());
    dVdZ.resize(valueDerivExpressions.size());
    dVdR1.resize(valueDerivExpressions.size());
    dVdR2.resize(valueDerivExpressions.size());
97
98
99
100
101
}

CpuCustomGBForce::~CpuCustomGBForce() {
}

102
void CpuCustomGBForce::setUseCutoff(float distance, const CpuNeighborList& neighbors) {
103
104
105
    cutoff = true;
    cutoffDistance = distance;
    neighborList = &neighbors;
106
  }
107
108
109
110
111
112
113
114
115
116
117
118
119
120

void CpuCustomGBForce::setPeriodic(RealVec& boxSize) {

    if (cutoff) {
        assert(boxSize[0] >= 2.0*cutoffDistance);
        assert(boxSize[1] >= 2.0*cutoffDistance);
        assert(boxSize[2] >= 2.0*cutoffDistance);
    }
    periodic = true;
    periodicBoxSize[0] = boxSize[0];
    periodicBoxSize[1] = boxSize[1];
    periodicBoxSize[2] = boxSize[2];
  }

121
122
123
124
125
void CpuCustomGBForce::calculateIxn(int numberOfAtoms, float* posq, RealOpenMM** atomParameters,
                                           const vector<set<int> >& exclusions, map<string, double>& globalParameters, float* forces,
                                           double* totalEnergy) {
    fvec4 boxSize(periodicBoxSize[0], periodicBoxSize[1], periodicBoxSize[2], 0);
    fvec4 invBoxSize((1/periodicBoxSize[0]), (1/periodicBoxSize[1]), (1/periodicBoxSize[2]), 0);
126
127
128
129
130
131
132
    for (map<string, double>::const_iterator iter = globalParameters.begin(); iter != globalParameters.end(); ++iter)
        expressionSet.setVariable(expressionSet.getVariableIndex(iter->first), iter->second);

    // First calculate the computed values.

    int numValues = valueTypes.size();
    for (int valueIndex = 0; valueIndex < numValues; valueIndex++) {
133
        if (valueTypes[valueIndex] == CustomGBForce::SingleParticle)
134
            calculateSingleParticleValue(valueIndex, numberOfAtoms, posq, values, atomParameters);
135
        else if (valueTypes[valueIndex] == CustomGBForce::ParticlePair)
136
            calculateParticlePairValue(valueIndex, numberOfAtoms, posq, atomParameters, values, exclusions, true, boxSize, invBoxSize);
137
        else
138
            calculateParticlePairValue(valueIndex, numberOfAtoms, posq, atomParameters, values, exclusions, false, boxSize, invBoxSize);
139
140
141
142
    }

    // Now calculate the energy and its derivatives.

143
144
145
    for (int i = 0; i < (int) dEdV.size(); i++)
        for (int j = 0; j < (int) dEdV[i].size(); j++)
            dEdV[i][j] = 0.0;
146
    for (int termIndex = 0; termIndex < (int) energyExpressions.size(); termIndex++) {
147
        if (energyTypes[termIndex] == CustomGBForce::SingleParticle)
148
            calculateSingleParticleEnergyTerm(termIndex, numberOfAtoms, posq, values, atomParameters, forces, totalEnergy, dEdV);
149
        else if (energyTypes[termIndex] == CustomGBForce::ParticlePair)
150
            calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, posq, atomParameters, values, exclusions, true, forces, totalEnergy, dEdV, boxSize, invBoxSize);
151
        else
152
            calculateParticlePairEnergyTerm(termIndex, numberOfAtoms, posq, atomParameters, values, exclusions, false, forces, totalEnergy, dEdV, boxSize, invBoxSize);
153
154
155
156
    }

    // Apply the chain rule to evaluate forces.

157
    calculateChainRuleForces(numberOfAtoms, posq, atomParameters, values, exclusions, forces, dEdV, boxSize, invBoxSize);
158
159
}

160
void CpuCustomGBForce::calculateSingleParticleValue(int index, int numAtoms, float* posq, vector<vector<float> >& values,
161
        RealOpenMM** atomParameters) {
162
163
    values[index].resize(numAtoms);
    for (int i = 0; i < numAtoms; i++) {
164
165
166
        expressionSet.setVariable(xindex, posq[4*i]);
        expressionSet.setVariable(yindex, posq[4*i+1]);
        expressionSet.setVariable(zindex, posq[4*i+2]);
167
168
169
170
        for (int j = 0; j < (int) paramNames.size(); j++)
            expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
        for (int j = 0; j < index; j++)
            expressionSet.setVariable(valueIndex[j], values[j][i]);
171
        values[index][i] = (float) valueExpressions[index].evaluate();
172
173
174
    }
}

175
176
void CpuCustomGBForce::calculateParticlePairValue(int index, int numAtoms, float* posq, RealOpenMM** atomParameters,
        vector<vector<float> >& values, const vector<set<int> >& exclusions, bool useExclusions, const fvec4& boxSize, const fvec4& invBoxSize) {
177
178
    values[index].resize(numAtoms);
    for (int i = 0; i < numAtoms; i++)
179
        values[index][i] = 0.0f;
180
181
182
    if (cutoff) {
        // Loop over all pairs in the neighbor list.

183
184
185
186
187
188
189
190
191
192
193
        for (int blockIndex = 0; blockIndex < neighborList->getNumBlocks(); blockIndex++) {
            const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
            const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
            const vector<char>& blockExclusions = neighborList->getBlockExclusions(blockIndex);
            for (int i = 0; i < (int) neighbors.size(); i++) {
                int first = neighbors[i];
                for (int k = 0; k < 4; k++) {
                    if ((blockExclusions[i] & (1<<k)) == 0) {
                        int second = blockAtom[k];
                        if (useExclusions && exclusions[first].find(second) != exclusions[first].end())
                            continue;
194
195
                        calculateOnePairValue(index, first, second, posq, atomParameters, values, boxSize, invBoxSize);
                        calculateOnePairValue(index, second, first, posq, atomParameters, values, boxSize, invBoxSize);
196
197
198
                    }
                }
            }
199
200
201
202
203
204
205
206
207
        }
    }
    else {
        // Perform an O(N^2) loop over all atom pairs.

        for (int i = 0; i < numAtoms; i++) {
            for (int j = i+1; j < numAtoms; j++) {
                if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
                    continue;
208
209
                calculateOnePairValue(index, i, j, posq, atomParameters, values, boxSize, invBoxSize);
                calculateOnePairValue(index, j, i, posq, atomParameters, values, boxSize, invBoxSize);
210
211
212
213
214
           }
        }
    }
}

215
216
217
218
219
220
221
222
void CpuCustomGBForce::calculateOnePairValue(int index, int atom1, int atom2, float* posq, RealOpenMM** atomParameters,
        vector<vector<float> >& values, const fvec4& boxSize, const fvec4& invBoxSize) {
    fvec4 deltaR;
    fvec4 pos1(posq+4*atom1);
    fvec4 pos2(posq+4*atom2);
    float r2;
    getDeltaR(pos2, pos1, deltaR, r2, periodic, boxSize, invBoxSize);
    float r = sqrtf(r2);
223
224
225
226
227
228
229
230
231
232
233
    if (cutoff && r >= cutoffDistance)
        return;
    for (int i = 0; i < (int) paramNames.size(); i++) {
        expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
        expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
    }
    expressionSet.setVariable(rindex, r);
    for (int i = 0; i < index; i++) {
        expressionSet.setVariable(particleValueIndex[i*2], values[i][atom1]);
        expressionSet.setVariable(particleValueIndex[i*2+1], values[i][atom2]);
    }
234
    values[index][atom1] += (float) valueExpressions[index].evaluate();
235
236
}

237
238
239
void CpuCustomGBForce::calculateSingleParticleEnergyTerm(int index, int numAtoms, float* posq, const vector<vector<float> >& values,
        RealOpenMM** atomParameters, float* forces, double* totalEnergy,
        vector<vector<float> >& dEdV) {
240
    for (int i = 0; i < numAtoms; i++) {
241
242
243
        expressionSet.setVariable(xindex, posq[4*i]);
        expressionSet.setVariable(yindex, posq[4*i+1]);
        expressionSet.setVariable(zindex, posq[4*i+2]);
244
245
246
247
248
        for (int j = 0; j < (int) paramNames.size(); j++)
            expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
        for (int j = 0; j < (int) valueNames.size(); j++)
            expressionSet.setVariable(valueIndex[j], values[j][i]);
        if (totalEnergy != NULL)
249
            *totalEnergy += (float) energyExpressions[index].evaluate();
250
        for (int j = 0; j < (int) valueNames.size(); j++)
251
252
253
254
            dEdV[j][i] += (float) energyDerivExpressions[index][j].evaluate();
        forces[4*i+0] -= (float) energyGradientExpressions[index][0].evaluate();
        forces[4*i+1] -= (float) energyGradientExpressions[index][1].evaluate();
        forces[4*i+2] -= (float) energyGradientExpressions[index][2].evaluate();
255
256
257
    }
}

258
259
260
void CpuCustomGBForce::calculateParticlePairEnergyTerm(int index, int numAtoms, float* posq, RealOpenMM** atomParameters,
        const vector<vector<float> >& values, const vector<set<int> >& exclusions, bool useExclusions,
        float* forces, double* totalEnergy, vector<vector<float> >& dEdV, const fvec4& boxSize, const fvec4& invBoxSize) {
261
262
263
    if (cutoff) {
        // Loop over all pairs in the neighbor list.

264
265
266
267
268
269
270
271
272
273
274
        for (int blockIndex = 0; blockIndex < neighborList->getNumBlocks(); blockIndex++) {
            const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
            const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
            const vector<char>& blockExclusions = neighborList->getBlockExclusions(blockIndex);
            for (int i = 0; i < (int) neighbors.size(); i++) {
                int first = neighbors[i];
                for (int k = 0; k < 4; k++) {
                    if ((blockExclusions[i] & (1<<k)) == 0) {
                        int second = blockAtom[k];
                        if (useExclusions && exclusions[first].find(second) != exclusions[first].end())
                            continue;
275
                        calculateOnePairEnergyTerm(index, first, second, posq, atomParameters, values, forces, totalEnergy, dEdV, boxSize, invBoxSize);
276
277
278
                    }
                }
            }
279
280
281
282
283
284
285
286
287
        }
    }
    else {
        // Perform an O(N^2) loop over all atom pairs.

        for (int i = 0; i < numAtoms; i++) {
            for (int j = i+1; j < numAtoms; j++) {
                if (useExclusions && exclusions[i].find(j) != exclusions[i].end())
                    continue;
288
                calculateOnePairEnergyTerm(index, i, j, posq, atomParameters, values, forces, totalEnergy, dEdV, boxSize, invBoxSize);
289
290
291
292
293
           }
        }
    }
}

294
295
296
void CpuCustomGBForce::calculateOnePairEnergyTerm(int index, int atom1, int atom2, float* posq, RealOpenMM** atomParameters,
        const vector<vector<float> >& values, float* forces, double* totalEnergy,
        vector<vector<float> >& dEdV, const fvec4& boxSize, const fvec4& invBoxSize) {
297
298
    // Compute the displacement.

299
300
301
302
303
304
    fvec4 deltaR;
    fvec4 pos1(posq+4*atom1);
    fvec4 pos2(posq+4*atom2);
    float r2;
    getDeltaR(pos2, pos1, deltaR, r2, periodic, boxSize, invBoxSize);
    float r = sqrtf(r2);
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
    if (cutoff && r >= cutoffDistance)
        return;

    // Record variables for evaluating expressions.

    for (int i = 0; i < (int) paramNames.size(); i++) {
        expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
        expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
    }
    expressionSet.setVariable(rindex, r);
    for (int i = 0; i < (int) valueNames.size(); i++) {
        expressionSet.setVariable(particleValueIndex[i*2], values[i][atom1]);
        expressionSet.setVariable(particleValueIndex[i*2+1], values[i][atom2]);
    }

    // Evaluate the energy and its derivatives.

    if (totalEnergy != NULL)
323
324
        *totalEnergy += (float) energyExpressions[index].evaluate();
    float dEdR = (float) energyDerivExpressions[index][0].evaluate();
325
    dEdR *= 1/r;
326
327
328
    fvec4 result = deltaR*dEdR;
    (fvec4(forces+4*atom1)-result).store(forces+4*atom1);
    (fvec4(forces+4*atom2)+result).store(forces+4*atom2);
329
    for (int i = 0; i < (int) valueNames.size(); i++) {
330
331
        dEdV[i][atom1] += (float) energyDerivExpressions[index][2*i+1].evaluate();
        dEdV[i][atom2] += (float) energyDerivExpressions[index][2*i+2].evaluate();
332
333
334
    }
}

335
336
337
void CpuCustomGBForce::calculateChainRuleForces(int numAtoms, float* posq, RealOpenMM** atomParameters,
        const vector<vector<float> >& values, const vector<set<int> >& exclusions, float* forces, vector<vector<float> >& dEdV,
        const fvec4& boxSize, const fvec4& invBoxSize) {
338
339
340
    if (cutoff) {
        // Loop over all pairs in the neighbor list.

341
342
343
344
345
346
347
348
349
350
        for (int blockIndex = 0; blockIndex < neighborList->getNumBlocks(); blockIndex++) {
            const int* blockAtom = &neighborList->getSortedAtoms()[4*blockIndex];
            const vector<int>& neighbors = neighborList->getBlockNeighbors(blockIndex);
            const vector<char>& blockExclusions = neighborList->getBlockExclusions(blockIndex);
            for (int i = 0; i < (int) neighbors.size(); i++) {
                int first = neighbors[i];
                for (int k = 0; k < 4; k++) {
                    if ((blockExclusions[i] & (1<<k)) == 0) {
                        int second = blockAtom[k];
                        bool isExcluded = (exclusions[first].find(second) != exclusions[first].end());
351
352
                        calculateOnePairChainRule(first, second, posq, atomParameters, values, forces, dEdV, isExcluded, boxSize, invBoxSize);
                        calculateOnePairChainRule(second, first, posq, atomParameters, values, forces, dEdV, isExcluded, boxSize, invBoxSize);
353
354
355
                    }
                }
            }
356
357
358
359
360
361
362
363
        }
    }
    else {
        // Perform an O(N^2) loop over all atom pairs.

        for (int i = 0; i < numAtoms; i++) {
            for (int j = i+1; j < numAtoms; j++) {
                bool isExcluded = (exclusions[i].find(j) != exclusions[i].end());
364
365
                calculateOnePairChainRule(i, j, posq, atomParameters, values, forces, dEdV, isExcluded, boxSize, invBoxSize);
                calculateOnePairChainRule(j, i, posq, atomParameters, values, forces, dEdV, isExcluded, boxSize, invBoxSize);
366
367
368
369
370
371
372
           }
        }
    }

    // Compute chain rule terms for computed values that depend explicitly on particle coordinates.

    for (int i = 0; i < numAtoms; i++) {
373
374
375
        expressionSet.setVariable(xindex, posq[4*i]);
        expressionSet.setVariable(yindex, posq[4*i+1]);
        expressionSet.setVariable(zindex, posq[4*i+2]);
376
377
378
379
        for (int j = 0; j < (int) paramNames.size(); j++)
            expressionSet.setVariable(paramIndex[j], atomParameters[i][j]);
        for (int j = 1; j < (int) valueNames.size(); j++) {
            expressionSet.setVariable(valueIndex[j-1], values[j-1][i]);
380
381
382
            dVdX[j] = 0.0;
            dVdY[j] = 0.0;
            dVdZ[j] = 0.0;
383
            for (int k = 1; k < j; k++) {
384
                float dVdV = (float) valueDerivExpressions[j][k].evaluate();
385
386
387
388
                dVdX[j] += dVdV*dVdX[k];
                dVdY[j] += dVdV*dVdY[k];
                dVdZ[j] += dVdV*dVdZ[k];
            }
389
390
391
392
393
394
            dVdX[j] += (float) valueGradientExpressions[j][0].evaluate();
            dVdY[j] += (float) valueGradientExpressions[j][1].evaluate();
            dVdZ[j] += (float) valueGradientExpressions[j][2].evaluate();
            forces[4*i+0] -= dEdV[j][i]*dVdX[j];
            forces[4*i+1] -= dEdV[j][i]*dVdY[j];
            forces[4*i+2] -= dEdV[j][i]*dVdZ[j];
395
396
397
398
        }
    }
}

399
400
void CpuCustomGBForce::calculateOnePairChainRule(int atom1, int atom2, float* posq, RealOpenMM** atomParameters,
        const vector<vector<float> >& values, float* forces, vector<vector<float> >& dEdV, bool isExcluded, const fvec4& boxSize, const fvec4& invBoxSize) {
401
402
    // Compute the displacement.

403
404
405
406
407
408
    fvec4 deltaR;
    fvec4 pos1(posq+4*atom1);
    fvec4 pos2(posq+4*atom2);
    float r2;
    getDeltaR(pos2, pos1, deltaR, r2, periodic, boxSize, invBoxSize);
    float r = sqrtf(r2);
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
    if (cutoff && r >= cutoffDistance)
        return;

    // Record variables for evaluating expressions.

    for (int i = 0; i < (int) paramNames.size(); i++) {
        expressionSet.setVariable(particleParamIndex[i*2], atomParameters[atom1][i]);
        expressionSet.setVariable(particleParamIndex[i*2+1], atomParameters[atom2][i]);
    }
    expressionSet.setVariable(rindex, r);
    expressionSet.setVariable(particleValueIndex[0], values[0][atom1]);
    expressionSet.setVariable(particleValueIndex[1], values[0][atom2]);

    // Evaluate the derivative of each parameter with respect to position and apply forces.

424
425
    float rinv = 1/r;
    deltaR *= rinv;
426
    if (!isExcluded || valueTypes[0] != CustomGBForce::ParticlePair) {
427
        dVdR1[0] = (float) valueDerivExpressions[0][0].evaluate();
428
        dVdR2[0] = -dVdR1[0];
429
430
        (fvec4(forces+4*atom1)-deltaR*(dEdV[0][atom1]*dVdR1[0])).store(forces+4*atom1);
        (fvec4(forces+4*atom2)-deltaR*(dEdV[0][atom1]*dVdR2[0])).store(forces+4*atom2);
431
432
433
434
435
436
    }
    for (int i = 0; i < (int) paramNames.size(); i++)
        expressionSet.setVariable(paramIndex[i], atomParameters[atom1][i]);
    expressionSet.setVariable(valueIndex[0], values[0][atom1]);
    for (int i = 1; i < (int) valueNames.size(); i++) {
        expressionSet.setVariable(valueIndex[i], values[i][atom1]);
437
438
439
        expressionSet.setVariable(xindex, posq[4*atom1]);
        expressionSet.setVariable(yindex, posq[4*atom1+1]);
        expressionSet.setVariable(zindex, posq[4*atom1+2]);
440
441
        dVdR1[i] = 0.0;
        dVdR2[i] = 0.0;
442
        for (int j = 0; j < i; j++) {
443
            float dVdV = (float) valueDerivExpressions[i][j].evaluate();
444
445
446
            dVdR1[i] += dVdV*dVdR1[j];
            dVdR2[i] += dVdV*dVdR2[j];
        }
447
448
449
450
451
452
453
454
455
456
        (fvec4(forces+4*atom1)-deltaR*(dEdV[i][atom1]*dVdR1[i])).store(forces+4*atom1);
        (fvec4(forces+4*atom2)-deltaR*(dEdV[i][atom1]*dVdR2[i])).store(forces+4*atom2);
    }
}

void CpuCustomGBForce::getDeltaR(const fvec4& posI, const fvec4& posJ, fvec4& deltaR, float& r2, bool periodic, const fvec4& boxSize, const fvec4& invBoxSize) const {
    deltaR = posJ-posI;
    if (periodic) {
        fvec4 base = round(deltaR*invBoxSize)*boxSize;
        deltaR = deltaR-base;
457
    }
458
    r2 = dot3(deltaR, deltaR);
459
}