TestReferenceAmoebaStretchBendForce.cpp 12.8 KB
Newer Older
Mark Friedrichs's avatar
Mark Friedrichs committed
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.      *
Mark Friedrichs's avatar
Mark Friedrichs committed
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
 * 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.                                     *
 * -------------------------------------------------------------------------- */

/**
 * This tests the Reference implementation of ReferenceAmoebaStretchBendForce.
 */

36
#include "openmm/internal/AssertionUtilities.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
37
38
//#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
39
#include "OpenMMAmoeba.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
40
41
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
42
#include "SimTKOpenMMRealType.h"
Mark Friedrichs's avatar
Mark Friedrichs committed
43
44
45
46
47
#include <iostream>
#include <vector>

using namespace OpenMM;

48
49
extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();

Mark Friedrichs's avatar
Mark Friedrichs committed
50
const double TOL = 1e-4;
51
const double DegreesToRadians = (M_PI/180.0);
Mark Friedrichs's avatar
Mark Friedrichs committed
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

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

   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

   --------------------------------------------------------------------------------------- */
     
67
static void crossProductVector3(double* vectorX, double* vectorY, double* vectorZ) {
Mark Friedrichs's avatar
Mark Friedrichs committed
68
69
70
71
72
73
74
75

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

76
static double dotVector3(double* vectorX, double* vectorY) {
Mark Friedrichs's avatar
Mark Friedrichs committed
77
78
79
80
81
    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
82
                                          std::vector<Vec3>& forces, double* energy) {
Mark Friedrichs's avatar
Mark Friedrichs committed
83
84

    int particle1, particle2, particle3;
85
    double abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend;
Mark Friedrichs's avatar
Mark Friedrichs committed
86

87
    amoebaStretchBendForce.getStretchBendParameters(bondIndex, particle1, particle2, particle3, abBondLength, cbBondLength, angleStretchBend, kStretchBend, k2StretchBend);
Mark Friedrichs's avatar
Mark Friedrichs committed
88
89
90
91
92
93
94
95
96
97
98
99
100
    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;
101
    for (int ii = 0; ii < 3; ii++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
102
103
104
105
106
107
         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];
    }
108
109
    double rAB   = sqrt(rAB2);
    double rCB   = sqrt(rCB2);
Mark Friedrichs's avatar
Mark Friedrichs committed
110

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

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

    // get positions and zero forces

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

    // calculates forces/energy

    *expectedEnergy = 0.0;
202
    for (int ii = 0; ii < amoebaStretchBendForce.getNumStretchBends(); ii++) {
peastman's avatar
peastman committed
203
        computeAmoebaStretchBendForce(ii, positions, amoebaStretchBendForce, expectedForces, expectedEnergy);
Mark Friedrichs's avatar
Mark Friedrichs committed
204
205
206
    }
}

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

    std::vector<Vec3> expectedForces;
    double expectedEnergy;
peastman's avatar
peastman committed
212
    computeAmoebaStretchBendForces(context, amoebaStretchBendForce, expectedForces, &expectedEnergy);
Mark Friedrichs's avatar
Mark Friedrichs committed
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);
Mark Friedrichs's avatar
Mark Friedrichs committed
218
    }
219
    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
Mark Friedrichs's avatar
Mark Friedrichs committed
220
221
}

peastman's avatar
peastman committed
222
void testOneStretchBend() {
Mark Friedrichs's avatar
Mark Friedrichs committed
223
224
225

    System system;
    int numberOfParticles = 3;
226
    for (int ii = 0; ii < numberOfParticles; ii++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
227
228
229
230
231
232
233
234
235
236
237
238
239
        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;

240
    amoebaStretchBendForce->addStretchBend(0, 1, 2, abLength, cbLength, angleStretchBend, kStretchBend, kStretchBend);
Mark Friedrichs's avatar
Mark Friedrichs committed
241
242

    system.addForce(amoebaStretchBendForce);
243
244
    ASSERT(!amoebaStretchBendForce->usesPeriodicBoundaryConditions());
    ASSERT(!system.usesPeriodicBoundaryConditions());
245
    Context context(system, integrator, Platform::getPlatformByName("Reference"));
Mark Friedrichs's avatar
Mark Friedrichs committed
246
247
248

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

249
250
251
    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);
Mark Friedrichs's avatar
Mark Friedrichs committed
252
253

    context.setPositions(positions);
peastman's avatar
peastman committed
254
    compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
255
256
257
    
    // Try changing the stretch-bend parameters and make sure it's still correct.
    
258
    amoebaStretchBendForce->setStretchBendParameters(0, 0, 1, 2, 1.1*abLength, 1.2*cbLength, 1.3*angleStretchBend, 1.4*kStretchBend, 1.4*kStretchBend);
259
260
261
    bool exceptionThrown = false;
    try {
        // This should throw an exception.
peastman's avatar
peastman committed
262
        compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
263
264
265
266
267
268
    }
    catch (std::exception ex) {
        exceptionThrown = true;
    }
    ASSERT(exceptionThrown);
    amoebaStretchBendForce->updateParametersInContext(context);
peastman's avatar
peastman committed
269
    compareWithExpectedForceAndEnergy(context, *amoebaStretchBendForce, TOL, "testOneStretchBend");
Mark Friedrichs's avatar
Mark Friedrichs committed
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
303
304
305
306
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("Reference"));
    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);
}

307
int main(int numberOfArguments, char* argv[]) {
Mark Friedrichs's avatar
Mark Friedrichs committed
308
309
310

    try {
        std::cout << "TestReferenceAmoebaStretchBendForce running test..." << std::endl;
311
        registerAmoebaReferenceKernelFactories();
peastman's avatar
peastman committed
312
        testOneStretchBend();
313
        testPeriodic();
Mark Friedrichs's avatar
Mark Friedrichs committed
314
315
316
317
318
319
    }
    catch(const std::exception& e) {
        std::cout << "exception: " << e.what() << std::endl;
        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
        return 1;
    }
320
321
    //std::cout << "PASS - Test succeeded." << std::endl;
    std::cout << "Done" << std::endl;
Mark Friedrichs's avatar
Mark Friedrichs committed
322
323
    return 0;
}