TestCudaAmoebaAngleForce.cpp 13.5 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
 * 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.                                     *
 * -------------------------------------------------------------------------- */

/**
33
 * This tests the CUDA implementation of CudaAmoebaAngleForce.
34
35
 */

peastman's avatar
peastman committed
36
37
38
#ifdef WIN32
  #define _USE_MATH_DEFINES // Needed to get M_PI
#endif
39
40
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
41
#include "openmm/CustomAngleForce.h"
42
43
44
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
45
#include "SimTKOpenMMRealType.h"
46
47
48
49
50
#include <iostream>
#include <vector>

using namespace OpenMM;

51
52
extern "C" void registerAmoebaCudaKernelFactories();

53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
const double TOL = 1e-5;

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

   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

   --------------------------------------------------------------------------------------- */
     
69
static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
70
71
72
73
74
75
76
77

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

78
79
80
static void getPrefactorsGivenAngleCosine(double cosine, double idealAngle, double quadraticK, double cubicK,
                                          double quarticK, double penticK, double sexticK,
                                          double* dEdR, double* energyTerm) {
81
82

    double angle;
83
    if (cosine >= 1.0) {
84
        angle = 0.0f;
85
    } else if (cosine <= -1.0) {
86
87
88
89
90
91
92
93
94
95
96
97
        angle = RADIAN*PI_M;
    } else {
        angle = RADIAN*acos(cosine);
    }

    double deltaIdeal         = angle - idealAngle;
    double deltaIdeal2        = deltaIdeal*deltaIdeal;
    double deltaIdeal3        = deltaIdeal*deltaIdeal2;
    double deltaIdeal4        = deltaIdeal2*deltaIdeal2;
 
    // deltaIdeal = r - r_0
 
98
99
100
101
102
    *dEdR        = (2.0                        +
                    3.0*cubicK*  deltaIdeal    +
                    4.0*quarticK*deltaIdeal2   +
                    5.0*penticK* deltaIdeal3   +
                    6.0*sexticK* deltaIdeal4    );
103
104
105
106
107
108
109
110
111
112
113
114
115
 
    *dEdR       *= RADIAN*quadraticK*deltaIdeal;
 

    *energyTerm  = 1.0f + cubicK* deltaIdeal    +
                          quarticK*deltaIdeal2   +
                          penticK* deltaIdeal3   +
                          sexticK* deltaIdeal4;
    *energyTerm *= quadraticK*deltaIdeal2;

    return;
}

116
static void computeAmoebaAngleForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaAngleForce& amoebaAngleForce,
117
                                             std::vector<Vec3>& forces, double* energy) {
118
119
120
121

    int particle1, particle2, particle3;
    double idealAngle;
    double quadraticK;
122
    amoebaAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, idealAngle, quadraticK);
123

124
125
126
127
    double cubicK         = amoebaAngleForce.getAmoebaGlobalAngleCubic();
    double quarticK       = amoebaAngleForce.getAmoebaGlobalAngleQuartic();
    double penticK        = amoebaAngleForce.getAmoebaGlobalAnglePentic();
    double sexticK        = amoebaAngleForce.getAmoebaGlobalAngleSextic();
128
129
130
131

    double deltaR[2][3];
    double r2_0 = 0.0;
    double r2_1 = 0.0;
132
    for (int ii = 0; ii < 3; ii++) {
133
134
135
136
137
138
139
140
141
142

           deltaR[0][ii]    = positions[particle1][ii] - positions[particle2][ii];
           r2_0            += deltaR[0][ii]*deltaR[0][ii];

           deltaR[1][ii]    = positions[particle3][ii] - positions[particle2][ii];
           r2_1            += deltaR[1][ii]*deltaR[1][ii];

    }

    double pVector[3];
143
144
145
    crossProductVector3(deltaR[0], deltaR[1], pVector);
    double rp      = sqrt(pVector[0]*pVector[0] + pVector[1]*pVector[1] + pVector[2]*pVector[2]);
    if (rp < 1.0e-06) {
146
147
148
149
150
151
152
       rp = 1.0e-06;
    }   
    double dot    = deltaR[0][0]*deltaR[1][0] + deltaR[0][1]*deltaR[1][1] + deltaR[0][2]*deltaR[1][2];
    double cosine = dot/sqrt(r2_0*r2_1);

    double dEdR;
    double energyTerm;
153
154
    getPrefactorsGivenAngleCosine(cosine, idealAngle, quadraticK, cubicK,
                                  quarticK, penticK, sexticK, &dEdR, &energyTerm);
155
156
157
158
159

    double termA  = -dEdR/(r2_0*rp);
    double termC  =  dEdR/(r2_1*rp);

    double deltaCrossP[3][3];
160
161
162
    crossProductVector3(deltaR[0], pVector, deltaCrossP[0]);
    crossProductVector3(deltaR[1], pVector, deltaCrossP[2]);
    for (int ii = 0; ii < 3; ii++) {
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
        deltaCrossP[0][ii] *= termA;
        deltaCrossP[2][ii] *= termC;
        deltaCrossP[1][ii]  = -1.0*(deltaCrossP[0][ii] + deltaCrossP[2][ii]);
    }

    forces[particle1][0]       += deltaCrossP[0][0];
    forces[particle1][1]       += deltaCrossP[0][1];
    forces[particle1][2]       += deltaCrossP[0][2];

    forces[particle2][0]       += deltaCrossP[1][0];
    forces[particle2][1]       += deltaCrossP[1][1];
    forces[particle2][2]       += deltaCrossP[1][2];

    forces[particle3][0]       += deltaCrossP[2][0];
    forces[particle3][1]       += deltaCrossP[2][1];
    forces[particle3][2]       += deltaCrossP[2][2];

    *energy                    += energyTerm;
}

183
184
static void computeAmoebaAngleForces(Context& context, AmoebaAngleForce& amoebaAngleForce,
                                             std::vector<Vec3>& expectedForces, double* expectedEnergy) {
185
186
187
188
189

    // get positions and zero forces

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

    // calculates forces/energy

    *expectedEnergy = 0.0;
199
200
    for (int ii = 0; ii < amoebaAngleForce.getNumAngles(); ii++) {
        computeAmoebaAngleForce(ii, positions, amoebaAngleForce, expectedForces, expectedEnergy);
201
202
203
204
205
206
    }

    return;

}

207
208
void compareWithExpectedForceAndEnergy(Context& context, AmoebaAngleForce& amoebaAngleForce,
                                       double tolerance, const std::string& idString) {
209
210
211

    std::vector<Vec3> expectedForces;
    double expectedEnergy;
212
    computeAmoebaAngleForces(context, amoebaAngleForce, expectedForces, &expectedEnergy);
213
214
215
216
   
    State state                      = context.getState(State::Forces | State::Energy);
    const std::vector<Vec3> forces   = state.getForces();

217
218
    for (unsigned int ii = 0; ii < forces.size(); ii++) {
        ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
219
    }
220
    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
221
222
}

peastman's avatar
peastman committed
223
void testOneAngle() {
224
225
226

    System system;
    int numberOfParticles = 3;
227
    for (int ii = 0; ii < numberOfParticles; ii++) {
228
229
230
231
232
        system.addParticle(1.0);
    }

    LangevinIntegrator integrator(0.0, 0.1, 0.01);

233
    AmoebaAngleForce* amoebaAngleForce = new AmoebaAngleForce();
234
235
236
237
238
239
240

    double angle      = 100.0;
    double quadraticK = 1.0;
    double cubicK     = 1.0e-01;
    double quarticK   = 1.0e-02;
    double penticK    = 1.0e-03;
    double sexticK    = 1.0e-04;
241
    amoebaAngleForce->addAngle(0, 1, 2, angle, quadraticK);
242

243
244
245
246
    amoebaAngleForce->setAmoebaGlobalAngleCubic(cubicK);
    amoebaAngleForce->setAmoebaGlobalAngleQuartic(quarticK);
    amoebaAngleForce->setAmoebaGlobalAnglePentic(penticK);
    amoebaAngleForce->setAmoebaGlobalAngleSextic(sexticK);
247

248
    system.addForce(amoebaAngleForce);
249
    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
250
251
252
253
254
255
256
257

    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);
258
    compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
259
260
261
262
263
264
265
    
    // Try changing the angle parameters and make sure it's still correct.
    
    amoebaAngleForce->setAngleParameters(0, 0, 1, 2, 1.1*angle, 1.4*quadraticK);
    bool exceptionThrown = false;
    try {
        // This should throw an exception.
266
        compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
267
268
269
270
271
272
    }
    catch (std::exception ex) {
        exceptionThrown = true;
    }
    ASSERT(exceptionThrown);
    amoebaAngleForce->updateParametersInContext(context);
273
    compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
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
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
void testPeriodic() {
    // Create a force that uses periodic boundary conditions, then compare to an identical custom force.
    
    System system;
    system.setDefaultPeriodicBoxVectors(Vec3(3, 0, 0), Vec3(0, 3, 0), Vec3(0, 0, 3));
    int numParticles = 3;
    for (int ii = 0; ii < numParticles; ii++)
        system.addParticle(1.0);
    LangevinIntegrator integrator(0.0, 0.1, 0.01);
    AmoebaAngleForce* amoebaAngleForce = new AmoebaAngleForce();
    double angle      = 100.0;
    double quadraticK = 1.0;
    double cubicK     = 1.0e-01;
    double quarticK   = 1.0e-02;
    double penticK    = 1.0e-03;
    double sexticK    = 1.0e-04;
    amoebaAngleForce->addAngle(0, 1, 2, angle, quadraticK);
    amoebaAngleForce->setAmoebaGlobalAngleCubic(cubicK);
    amoebaAngleForce->setAmoebaGlobalAngleQuartic(quarticK);
    amoebaAngleForce->setAmoebaGlobalAnglePentic(penticK);
    amoebaAngleForce->setAmoebaGlobalAngleSextic(sexticK);
    amoebaAngleForce->setUsesPeriodicBoundaryConditions(true);
    system.addForce(amoebaAngleForce);
    CustomAngleForce* customForce = new CustomAngleForce("k2*delta^2 + k3*delta^3 + k4*delta^4 + k5*delta^5 + k6*delta^6; delta=theta-theta0");
    customForce->addGlobalParameter("theta0", angle*M_PI/180);
    customForce->addGlobalParameter("k2", quadraticK*pow(180/M_PI, 2.0));
    customForce->addGlobalParameter("k3", cubicK*pow(180/M_PI, 3.0));
    customForce->addGlobalParameter("k4", quarticK*pow(180/M_PI, 4.0));
    customForce->addGlobalParameter("k5", penticK*pow(180/M_PI, 5.0));
    customForce->addGlobalParameter("k6", sexticK*pow(180/M_PI, 6.0));
    customForce->addAngle(0, 1, 2);
    customForce->setUsesPeriodicBoundaryConditions(true);
    customForce->setForceGroup(1);
    system.addForce(customForce);
    Context context(system, integrator, Platform::getPlatformByName("CUDA"));

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

    positions[0] = Vec3(0, 1, 0);
    positions[1] = Vec3(0, 0, 0);
    positions[2] = Vec3(0, 0, 2);

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

326
int main(int argc, char* argv[]) {
327
    try {
328
        std::cout << "TestCudaAmoebaAngleForce running test..." << std::endl;
329
        registerAmoebaCudaKernelFactories();
330
        if (argc > 1)
331
            Platform::getPlatformByName("CUDA").setPropertyDefaultValue("Precision", std::string(argv[1]));
peastman's avatar
peastman committed
332
        testOneAngle();
333
        testPeriodic();
334
335
336
337
338
339
340
341
342

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