TestReferenceAmoebaBondForce.cpp 9.8 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
/* -------------------------------------------------------------------------- *
 *                                   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.               *
 *                                                                            *
 * Portions copyright (c) 2008 Stanford University and the Authors.           *
 * 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 Reference implementation of AmoebaBondForce.
34
35
 */

36
#include "openmm/internal/AssertionUtilities.h"
37
#include "openmm/Context.h"
38
#include "OpenMMAmoeba.h"
39
40
41
42
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include <iostream>
#include <vector>
43
#include <math.h>
44
45
46

using namespace OpenMM;

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

49
50
const double TOL = 1e-5;

51
static void computeAmoebaBondForce(int bondIndex,  std::vector<Vec3>& positions, AmoebaBondForce& AmoebaBondForce,
52
                                           std::vector<Vec3>& forces, double* energy) {
53
54
55
56

    int particle1, particle2;
    double bondLength;
    double quadraticK;
57
58
    double cubicK    = AmoebaBondForce.getAmoebaGlobalBondCubic();
    double quarticK  = AmoebaBondForce.getAmoebaGlobalBondQuartic();
59
    AmoebaBondForce.getBondParameters(bondIndex, particle1, particle2,  bondLength,  quadraticK);
60
61
62

    double deltaR[3];
    double r2 = 0.0;
63
    for (int ii = 0; ii < 3; ii++) {
64
65
66
           deltaR[ii]    = positions[particle2][ii] - positions[particle1][ii];
           r2           += deltaR[ii]*deltaR[ii];
    }
67
    double r                   = sqrt(r2);
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

    double bondDelta           = (r - bondLength);
    double bondDelta2          = bondDelta*bondDelta;
    double dEdR                = 1.0 + 1.5*cubicK*bondDelta + 2.0*quarticK*bondDelta2;

           dEdR               *= (r > 0.0) ? (2.0*quadraticK*bondDelta)/r : 0.0;

   forces[particle1][0]       += dEdR*deltaR[0];
   forces[particle1][1]       += dEdR*deltaR[1];
   forces[particle1][2]       += dEdR*deltaR[2];

   forces[particle2][0]       -= dEdR*deltaR[0];
   forces[particle2][1]       -= dEdR*deltaR[1];
   forces[particle2][2]       -= dEdR*deltaR[2];

   *energy                    += (1.0f + cubicK*bondDelta + quarticK*bondDelta2)*quadraticK*bondDelta2;

}

87
88
static void computeAmoebaBondForces(Context& context, AmoebaBondForce& AmoebaBondForce,
                                             std::vector<Vec3>& expectedForces, double* expectedEnergy, FILE* log) {
89
90
91
92
93

    // get positions and zero forces

    State state = context.getState(State::Positions);
    std::vector<Vec3> positions = state.getPositions();
94
    expectedForces.resize(positions.size());
95
    
96
    for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
97
98
99
100
101
102
        expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
    }

    // calculates forces/energy

    *expectedEnergy = 0.0;
103
104
    for (int ii = 0; ii < AmoebaBondForce.getNumBonds(); ii++) {
        computeAmoebaBondForce(ii, positions, AmoebaBondForce, expectedForces, expectedEnergy);
105
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
106
#ifdef AMOEBA_DEBUG
107
108
109
110
    if (log) {
        (void) fprintf(log, "computeAmoebaBondForces: expected energy=%15.7e\n", *expectedEnergy);
        for (unsigned int ii = 0; ii < positions.size(); ii++) {
            (void) fprintf(log, "%6u [%15.7e %15.7e %15.7e]\n", ii, expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2]);
111
        }
112
        (void) fflush(log);
113
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
114
#endif
115
116
117
118
    return;

}

119
void compareWithExpectedForceAndEnergy(Context& context, AmoebaBondForce& AmoebaBondForce, double tolerance, const std::string& idString, FILE* log) {
120
121
122

    std::vector<Vec3> expectedForces;
    double expectedEnergy;
123
    computeAmoebaBondForces(context, AmoebaBondForce, expectedForces, &expectedEnergy, NULL);
124
125
126
   
    State state                      = context.getState(State::Forces | State::Energy);
    const std::vector<Vec3> forces   = state.getForces();
Mark Friedrichs's avatar
Mark Friedrichs committed
127
#ifdef AMOEBA_DEBUG
128
129
130
131
132
    if (log) {
        (void) fprintf(log, "computeAmoebaBondForces: expected energy=%15.7e %15.7e\n", expectedEnergy, state.getPotentialEnergy());
        for (unsigned int ii = 0; ii < forces.size(); ii++) {
            (void) fprintf(log, "%6u [%15.7e %15.7e %15.7e]   [%15.7e %15.7e %15.7e]\n", ii,
                            expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2]);
133
        }
134
        (void) fflush(log);
135
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
136
#endif
137

138
139
    for (unsigned int ii = 0; ii < forces.size(); ii++) {
        ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
140
    }
141
    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
142
143
}

144
void testOneBond(FILE* log) {
145
146
147
148
149
150
151
152

    System system;

    system.addParticle(1.0);
    system.addParticle(1.0);

    LangevinIntegrator integrator(0.0, 0.1, 0.01);

153
    AmoebaBondForce* amoebaBondForce = new AmoebaBondForce();
154
155
156
157
158

    double bondLength = 1.5;
    double quadraticK = 1.0;
    double cubicK     = 2.0;
    double quarticicK = 3.0;
159
160
    amoebaBondForce->setAmoebaGlobalBondCubic(cubicK);
    amoebaBondForce->setAmoebaGlobalBondQuartic(quarticicK);
161
    amoebaBondForce->addBond(0, 1, bondLength, quadraticK);
162

163
    system.addForce(amoebaBondForce);
164
    Context context(system, integrator, Platform::getPlatformByName("Reference"));
165
166
167
168
169
170
    std::vector<Vec3> positions(2);

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

    context.setPositions(positions);
171
    compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testOneBond", log);
172
173
}

174
void testTwoBond(FILE* log) {
175
176
177
178
179
180
181
182
183

    System system;

    system.addParticle(1.0);
    system.addParticle(1.0);
    system.addParticle(1.0);

    LangevinIntegrator integrator(0.0, 0.1, 0.01);

184
    AmoebaBondForce* amoebaBondForce = new AmoebaBondForce();
185
186
187
188
189

    double bondLength = 1.5;
    double quadraticK = 1.0;
    double cubicK     = 2.0;
    double quarticicK = 3.0;
190
191
    amoebaBondForce->setAmoebaGlobalBondCubic(cubicK);
    amoebaBondForce->setAmoebaGlobalBondQuartic(quarticicK);
192
193
    amoebaBondForce->addBond(0, 1, bondLength, quadraticK);
    amoebaBondForce->addBond(1, 2, bondLength, quadraticK);
194

195
    system.addForce(amoebaBondForce);
196
197
    ASSERT(!amoebaBondForce->usesPeriodicBoundaryConditions());
    ASSERT(!system.usesPeriodicBoundaryConditions());
198
    Context context(system, integrator, Platform::getPlatformByName("Reference"));
199
200
201
202
203
204
205
    std::vector<Vec3> positions(3);

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

    context.setPositions(positions);
206
    compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond", log);
207
208
209
210
211
212
213
214
    
    // Try changing the bond parameters and make sure it's still correct.
    
    amoebaBondForce->setBondParameters(0, 0, 1, 1.1*bondLength, 1.4*quadraticK);
    amoebaBondForce->setBondParameters(1, 1, 2, 1.2*bondLength, 0.9*quadraticK);
    bool exceptionThrown = false;
    try {
        // This should throw an exception.
215
        compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond", log);
216
217
218
219
220
221
    }
    catch (std::exception ex) {
        exceptionThrown = true;
    }
    ASSERT(exceptionThrown);
    amoebaBondForce->updateParametersInContext(context);
222
    compareWithExpectedForceAndEnergy(context, *amoebaBondForce, TOL, "testTwoBond", log);
223
224
}

225
int main(int numberOfArguments, char* argv[]) {
226
227

    try {
228
        std::cout << "TestReferenceAmoebaBondForce running test..." << std::endl;
229
        registerAmoebaReferenceKernelFactories();
Mark Friedrichs's avatar
Mark Friedrichs committed
230
231
        FILE* log = NULL;
        //FILE* log = stderr;
232

233
234
        //testOneBond(log);
        testTwoBond(log);
Mark Friedrichs's avatar
Mark Friedrichs committed
235
#ifdef AMOEBA_DEBUG
236
237
        if (log && log != stderr)
            (void) fclose(log);
Mark Friedrichs's avatar
Mark Friedrichs committed
238
#endif
239
240
241
242
243
244
245

    }
    catch(const std::exception& e) {
        std::cout << "exception: " << e.what() << std::endl;
        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
        return 1;
    }
246
247
    //std::cout << "PASS - Test succeeded." << std::endl;
    std::cout << "Done" << std::endl;
248
249
    return 0;
}