TestCudaAmoebaStretchBendForce.cpp 12.7 KB
Newer Older
1
2
3
4
5
6
7
8
/* -------------------------------------------------------------------------- *
 *                                   OpenMMAmoeba                             *
 * -------------------------------------------------------------------------- *
 * 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
32
33
34
35
36
37
38
39
40
 * Authors: Mark Friedrichs                                                   *
 * 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.                                     *
 * -------------------------------------------------------------------------- */

/**
 * This tests the CUDA implementation of CudaAmoebaStretchBendForce.
 */

#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
41
#include "SimTKOpenMMRealType.h"
42
43
44
45
46
#include <iostream>
#include <vector>

using namespace OpenMM;

47
48
extern "C" void registerAmoebaCudaKernelFactories();

49
const double TOL = 1e-4;
50
const double DegreesToRadians = PI_M/180.0;
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

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

   Compute cross product of two 3-vectors and place in 3rd vector

   vectorZ = vectorX x vectorY

   @param vectorX             x-vector
   @param vectorY             y-vector
   @param vectorZ             z-vector

   @return vector is vectorZ

   --------------------------------------------------------------------------------------- */
     
66
static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
67
68
69
70
71
72
73
74

    vectorZ[0]  = vectorX[1]*vectorY[2] - vectorX[2]*vectorY[1];
    vectorZ[1]  = vectorX[2]*vectorY[0] - vectorX[0]*vectorY[2];
    vectorZ[2]  = vectorX[0]*vectorY[1] - vectorX[1]*vectorY[0];

    return;
}

75
static double dotVector3(double* vectorX, double* vectorY) {
76
77
78
79
80
    return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
}


static void computeAmoebaStretchBendForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaStretchBendForce& amoebaStretchBendForce,
peastman's avatar
peastman committed
81
                                          std::vector<Vec3>& forces, double* energy) {
82
83

    int particle1, particle2, particle3;
84
    double abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend;
85

86
    amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend);
87
88
89
90
91
92
93
94
95
96
97
98
    angleStretchBend *= RADIAN;
    enum { A, B, C, LastAtomIndex };
    enum { AB, CB, CBxAB, ABxP, CBxP, LastDeltaIndex };
 
    // ---------------------------------------------------------------------------------------
 
    // get deltaR between various combinations of the 3 atoms
    // and various intermediate terms
 
    double deltaR[LastDeltaIndex][3];
    double rAB2 = 0.0;
    double rCB2 = 0.0;
99
    for (int ii = 0; ii < 3; ii++) {
100
101
102
103
104
105
         deltaR[AB][ii]  = positions[particle1][ii] - positions[particle2][ii];
         rAB2           += deltaR[AB][ii]*deltaR[AB][ii];

         deltaR[CB][ii]  = positions[particle3][ii] - positions[particle2][ii];
         rCB2           += deltaR[CB][ii]*deltaR[CB][ii];
    }
106
107
    double rAB   = sqrt(rAB2);
    double rCB   = sqrt(rCB2);
108

109
110
111
    crossProductVector3(deltaR[CB], deltaR[AB], deltaR[CBxAB]);
    double  rP   = dotVector3(deltaR[CBxAB], deltaR[CBxAB]);
            rP   = sqrt(rP);
112
 
113
    if (rP <= 0.0) {
114
115
       return;
    }
116
    double dot    = dotVector3(deltaR[CB], deltaR[AB]);
117
118
119
    double cosine = dot/(rAB*rCB);
 
    double angle;
120
    if (cosine >= 1.0) {
121
       angle = 0.0;
122
    }
123
    else if (cosine <= -1.0) {
124
       angle = PI_M;
125
126
    }
    else {
127
128
129
130
131
132
133
134
       angle = RADIAN*acos(cosine);
    }
 
    double termA = -RADIAN/(rAB2*rP);
    double termC =  RADIAN/(rCB2*rP);
 
    // P = CBxAB
 
135
136
137
    crossProductVector3(deltaR[AB], deltaR[CBxAB], deltaR[ABxP]);
    crossProductVector3(deltaR[CB], deltaR[CBxAB], deltaR[CBxP]);
    for (int ii = 0; ii < 3; ii++) {
138
139
140
141
       deltaR[ABxP][ii] *= termA;
       deltaR[CBxP][ii] *= termC;
    }
 
142
143
    double dr1   = rAB - abBondLength;
    double dr2   = rCB - cbBondLength;
144
145
146
147
 
    termA        = 1.0/rAB;
    termC        = 1.0/rCB;
 
148
    double drkk = dr1 * kStretchBend + dr2 * k2StretchBend;
149
150
151
152
153
154
 
    // ---------------------------------------------------------------------------------------
 
    // forces
 
    // calculate forces for atoms a, b, c
155
    // the force for b is then -(a + c)
156
157
158
 
    double subForce[LastAtomIndex][3];
    double dt = angle - angleStretchBend;
159
    for (int jj = 0; jj < 3; jj++) {
160
161
        subForce[A][jj] = kStretchBend*dt*termA*deltaR[AB][jj] + drkk*deltaR[ABxP][jj];
        subForce[C][jj] = k2StretchBend*dt*termC*deltaR[CB][jj] + drkk*deltaR[CBxP][jj];
162
        subForce[B][jj] = -(subForce[A][jj] + subForce[C][jj]);
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
    }
 
    // ---------------------------------------------------------------------------------------
 
    // accumulate forces and energy
 
    forces[particle1][0]       -= subForce[0][0];
    forces[particle1][1]       -= subForce[0][1];
    forces[particle1][2]       -= subForce[0][2];

    forces[particle2][0]       -= subForce[1][0];
    forces[particle2][1]       -= subForce[1][1];
    forces[particle2][2]       -= subForce[1][2];

    forces[particle3][0]       -= subForce[2][0];
    forces[particle3][1]       -= subForce[2][1];
    forces[particle3][2]       -= subForce[2][2];

181
    *energy                    += dt*drkk;
182
183
}
 
184
static void computeAmoebaStretchBendForces(Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
peastman's avatar
peastman committed
185
                                          std::vector<Vec3>& expectedForces, double* expectedEnergy) {
186
187
188
189
190

    // get positions and zero forces

    State state                 = context.getState(State::Positions);
    std::vector<Vec3> positions = state.getPositions();
191
    expectedForces.resize(positions.size());
192
    
193
    for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
194
195
196
197
198
199
        expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
    }

    // calculates forces/energy

    *expectedEnergy = 0.0;
200
201
    for (int ii = 0; ii < amoebaStretchBendForce.getNumStretchBends(); ii++) {
        computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy);
202
203
204
    }
}

205
void compareWithExpectedForceAndEnergy(Context& context, AmoebaStretchBendForce& amoebaStretchBendForce,
peastman's avatar
peastman committed
206
                                        double tolerance, const std::string& idString) {
207
208
209

    std::vector<Vec3> expectedForces;
    double expectedEnergy;
210
    computeAmoebaStretchBendForces(context, amoebaStretchBendForce, expectedForces, &expectedEnergy);
211
212
213
   
    State state                      = context.getState(State::Forces | State::Energy);
    const std::vector<Vec3> forces   = state.getForces();
214
215
    for (unsigned int ii = 0; ii < forces.size(); ii++) {
        ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
216
    }
217
    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
218
219
}

peastman's avatar
peastman committed
220
void testOneStretchBend() {
221
222
223

    System system;
    int numberOfParticles = 3;
224
    for (int ii = 0; ii < numberOfParticles; ii++) {
225
226
227
228
229
230
231
232
233
234
235
236
237
        system.addParticle(1.0);
    }

    LangevinIntegrator integrator(0.0, 0.1, 0.01);

    AmoebaStretchBendForce* amoebaStretchBendForce = new AmoebaStretchBendForce();

    double abLength         = 0.144800000E+01;
    double cbLength         = 0.101500000E+01;
    double angleStretchBend = 0.108500000E+03*DegreesToRadians;
    //double kStretchBend     = 0.750491578E-01;
    double kStretchBend     = 1.0;

238
    amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend);
239
240

    system.addForce(amoebaStretchBendForce);
241
    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
242
243
244

    std::vector<Vec3> positions(numberOfParticles);

245
246
247
    positions[0] = Vec3(0.262660000E+02,  0.254130000E+02,  0.284200000E+01);
    positions[1] = Vec3(0.273400000E+02,  0.244300000E+02,  0.261400000E+01);
    positions[2] = Vec3(0.269573220E+02,  0.236108860E+02,  0.216376800E+01);
248
249

    context.setPositions(positions);
250
    compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
251
252
253
    
    // Try changing the stretch-bend parameters and make sure it's still correct.
    
254
    amoebaStretchBendForce->setStretchBendParameters(0, 0, 1, 2, 1.1*abLength, 1.2*cbLength, 1.3*angleStretchBend, 1.4*kStretchBend, 1.4*kStretchBend);
255
256
257
    bool exceptionThrown = false;
    try {
        // This should throw an exception.
258
        compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
259
260
261
262
263
264
    }
    catch (std::exception ex) {
        exceptionThrown = true;
    }
    ASSERT(exceptionThrown);
    amoebaStretchBendForce->updateParametersInContext(context);
265
    compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
266
267
}

268
269
270
271
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
void testPeriodic() {
    // Create a force that uses periodic boundary conditions.

    System system;
    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
    int numberOfParticles = 3;
    for (int ii = 0; ii < numberOfParticles; ii++)
        system.addParticle(1.0);
    LangevinIntegrator integrator(0.0, 0.1, 0.01);
    AmoebaStretchBendForce* amoebaStretchBendForce = new AmoebaStretchBendForce();
    double abLength         = 0.144800000E+01;
    double cbLength         = 0.101500000E+01;
    double angleStretchBend = 0.108500000E+03*DegreesToRadians;
    double kStretchBend     = 1.0;
    amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend);
    amoebaStretchBendForce->setUsesPeriodicBoundaryConditions(true);
    system.addForce(amoebaStretchBendForce);
    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
    std::vector<Vec3> positions(numberOfParticles);
    positions[0] = Vec3(0, 1, 0);
    positions[1] = Vec3(0, 0, 0);
    positions[2] = Vec3(0, 0, 1);
    context.setPositions(positions);
    State s1 = context.getState(State::Forces | State::Energy);
    
    // Move one atom to a position that should give identical results.

    positions[2] = Vec3(0, 0, -2);
    context.setPositions(positions);
    State s2 = context.getState(State::Forces | State::Energy);
    ASSERT_EQUAL_TOL(s1.getPotentialEnergy(), s2.getPotentialEnergy(), 1e-5);
    for (int i = 0; i < numberOfParticles; i++)
        ASSERT_EQUAL_VEC(s1.getForces()[i], s2.getForces()[i], 1e-5);
}

303
int main(int argc, char* argv[]) {
304
305
306
    try {
        std::cout << "TestCudaAmoebaStretchBendForce running test..." << std::endl;
        registerAmoebaCudaKernelFactories();
307
        if (argc > 1)
308
            Platform::getPlatformByName("CUDA").setPropertyDefaultValue("Precision", std::string(argv[1]));
peastman's avatar
peastman committed
309
        testOneStretchBend();
310
        testPeriodic();
311
312
313
314
315
316
317
318
319
320
321
    } catch(const std::exception& e) {
        std::cout << "exception: " << e.what() << std::endl;
        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
        return 1;
    }
    std::cout << "Done" << std::endl;
    return 0;
}