"devtools/packaging/scripts/vscode:/vscode.git/clone" did not exist on "bf6d95c2abc400d1af2e508b26aeb913a1fbc00b"
TestReferenceAmoebaAngleForce.cpp 13.6 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
32
 * 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.                                     *
 * -------------------------------------------------------------------------- */

/**
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
#include "openmm/internal/AssertionUtilities.h"
40
#include "openmm/Context.h"
41
#include "openmm/CustomAngleForce.h"
42
#include "OpenMMAmoeba.h"
43
44
#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" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();

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
static void getPrefactorsGivenAngleCosine(double cosine, double idealAngle, double quadraticK, double cubicK,
                                          double quarticK, double penticK, double sexticK,
peastman's avatar
peastman committed
80
                                          double* dEdR, double* energyTerm) {
81
82

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

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

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

    return;
}

118
static void computeAmoebaAngleForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaAngleForce& AmoebaAngleForce,
peastman's avatar
peastman committed
119
                                             std::vector<Vec3>& forces, double* energy) {
120
121
122
123

    int particle1, particle2, particle3;
    double idealAngle;
    double quadraticK;
124
    AmoebaAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, idealAngle, quadraticK);
125

126
127
128
129
    double cubicK         = AmoebaAngleForce.getAmoebaGlobalAngleCubic();
    double quarticK       = AmoebaAngleForce.getAmoebaGlobalAngleQuartic();
    double penticK        = AmoebaAngleForce.getAmoebaGlobalAnglePentic();
    double sexticK        = AmoebaAngleForce.getAmoebaGlobalAngleSextic();
130
131
132
133

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

           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];
145
146
147
    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) {
148
       rp = 1.0e-06;
Jason Swails's avatar
More  
Jason Swails committed
149
    }
150
151
152
153
154
    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;
155
    getPrefactorsGivenAngleCosine(cosine, idealAngle, quadraticK, cubicK,
peastman's avatar
peastman committed
156
                                  quarticK, penticK, sexticK, &dEdR, &energyTerm);
157
158
159
160
161

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

    double deltaCrossP[3][3];
162
163
164
    crossProductVector3(deltaR[0], pVector, deltaCrossP[0]);
    crossProductVector3(deltaR[1], pVector, deltaCrossP[2]);
    for (int ii = 0; ii < 3; ii++) {
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
        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;
}

185
static void computeAmoebaAngleForces(Context& context, AmoebaAngleForce& AmoebaAngleForce,
peastman's avatar
peastman committed
186
                                             std::vector<Vec3>& expectedForces, double* expectedEnergy) {
187
188
189
190
191

    // get positions and zero forces

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

    // calculates forces/energy

    *expectedEnergy = 0.0;
201
    for (int ii = 0; ii < AmoebaAngleForce.getNumAngles(); ii++) {
peastman's avatar
peastman committed
202
        computeAmoebaAngleForce(ii, positions, AmoebaAngleForce, expectedForces, expectedEnergy);
203
204
205
    }
}

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

    std::vector<Vec3> expectedForces;
    double expectedEnergy;
peastman's avatar
peastman committed
211
    computeAmoebaAngleForces(context, AmoebaAngleForce, expectedForces, &expectedEnergy);
212
213
214
215
   
    State state                      = context.getState(State::Forces | State::Energy);
    const std::vector<Vec3> forces   = state.getForces();

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

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

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

    LangevinIntegrator integrator(0.0, 0.1, 0.01);

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

    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;
240
    amoebaAngleForce->addAngle(0, 1, 2, angle, quadraticK);
241

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

247
    system.addForce(amoebaAngleForce);
248
249
    ASSERT(!amoebaAngleForce->usesPeriodicBoundaryConditions());
    ASSERT(!system.usesPeriodicBoundaryConditions());
250
    Context context(system, integrator, Platform::getPlatformByName("Reference"));
251
252
253
254
255
256
257
258

    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);
peastman's avatar
peastman committed
259
    compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
260
261
262
263
264
265
266
    
    // 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.
peastman's avatar
peastman committed
267
        compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
268
269
270
271
272
273
    }
    catch (std::exception ex) {
        exceptionThrown = true;
    }
    ASSERT(exceptionThrown);
    amoebaAngleForce->updateParametersInContext(context);
peastman's avatar
peastman committed
274
    compareWithExpectedForceAndEnergy(context, *amoebaAngleForce, TOL, "testOneAngle");
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
326
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("Reference"));

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

327
int main(int numberOfArguments, char* argv[]) {
328
329

    try {
330
        std::cout << "TestReferenceAmoebaAngleForce running test..." << std::endl;
331
        registerAmoebaReferenceKernelFactories();
peastman's avatar
peastman committed
332
        testOneAngle();
333
        testPeriodic();
334
335
336
337
338
339
    }
    catch(const std::exception& e) {
        std::cout << "exception: " << e.what() << std::endl;
        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
        return 1;
    }
340
341
    //std::cout << "PASS - Test succeeded." << std::endl;
    std::cout << "Done" << std::endl;
342
343
    return 0;
}