TestReferenceAmoebaInPlaneAngleForce.cpp 16 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
 * 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 ReferenceAmoebaInPlaneAngleForce.
Mark Friedrichs's avatar
Mark Friedrichs committed
34
35
 */

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

using namespace OpenMM;

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

Mark Friedrichs's avatar
Mark Friedrichs committed
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
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

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

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

74
static double dotVector3(double* vectorX, double* vectorY) {
Mark Friedrichs's avatar
Mark Friedrichs committed
75
76
77
    return vectorX[0]*vectorY[0] + vectorX[1]*vectorY[1] + vectorX[2]*vectorY[2];
}

78
79
static void getPrefactorsGivenInPlaneAngleCosine(double cosine, double idealInPlaneAngle, double quadraticK, double cubicK,
                                                 double quarticK, double penticK, double sexticK,
peastman's avatar
peastman committed
80
                                                 double* dEdR, double* energyTerm) {
Mark Friedrichs's avatar
Mark Friedrichs committed
81
82

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

    double deltaIdeal         = angle - idealInPlaneAngle;
    double deltaIdeal2        = deltaIdeal*deltaIdeal;
    double deltaIdeal3        = deltaIdeal*deltaIdeal2;
    double deltaIdeal4        = deltaIdeal2*deltaIdeal2;
 
    // deltaIdeal = r - r_0
 
100
    *dEdR        =  (2.0                        +
Mark Friedrichs's avatar
Mark Friedrichs committed
101
102
103
                     3.0*cubicK*  deltaIdeal    +
                     4.0*quarticK*deltaIdeal2   +
                     5.0*penticK* deltaIdeal3   +
104
                     6.0*sexticK* deltaIdeal4   );
Mark Friedrichs's avatar
Mark Friedrichs committed
105
106
107
108
 
    *dEdR       *= RADIAN*quadraticK*deltaIdeal;
 

109
110
111
    *energyTerm  = 1.0f + cubicK* deltaIdeal   +
                          quarticK*deltaIdeal2 +
                          penticK* deltaIdeal3 +
Mark Friedrichs's avatar
Mark Friedrichs committed
112
113
114
115
116
117
                          sexticK* deltaIdeal4;
    *energyTerm *= quadraticK*deltaIdeal2;

    return;
}

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

    int particle1, particle2, particle3, particle4;
    double idealInPlaneAngle;
    double quadraticK;
124
    AmoebaInPlaneAngleForce.getAngleParameters(bondIndex, particle1, particle2, particle3, particle4, idealInPlaneAngle, quadraticK);
Mark Friedrichs's avatar
Mark Friedrichs committed
125

126
127
128
129
    double cubicK         = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleCubic();
    double quarticK       = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleQuartic();
    double penticK        = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAnglePentic();
    double sexticK        = AmoebaInPlaneAngleForce.getAmoebaGlobalInPlaneAngleSextic();
Mark Friedrichs's avatar
Mark Friedrichs committed
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149

    // T   = AD x CD
    // P   = B + T*delta
    // AP  = A - P
    // CP  = A - P
    // M   = CP x AP

    enum { AD, BD, CD, T, AP, P, CP, M, APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex };
 
    // AD   0
    // BD   1
    // CD   2 
    //  T   3
    // AP   4
    //  P   5
    // CP   6
    // M    7
    // APxM, CPxM, ADxBD, BDxCD, TxCD, ADxT, dBxAD, CDxdB, LastDeltaAtomIndex

    double deltaR[LastDeltaAtomIndex][3];
150
    for (int ii = 0; ii < 3; ii++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
151
152
153
154
        deltaR[AD][ii] = positions[particle1][ii] - positions[particle4][ii];
        deltaR[BD][ii] = positions[particle2][ii] - positions[particle4][ii];
        deltaR[CD][ii] = positions[particle3][ii] - positions[particle4][ii];
    }
155
    crossProductVector3(deltaR[AD], deltaR[CD], deltaR[T]);
Mark Friedrichs's avatar
Mark Friedrichs committed
156
 
157
158
    double rT2     = dotVector3(deltaR[T], deltaR[T]);
    double delta   = dotVector3(deltaR[T], deltaR[BD]);
Mark Friedrichs's avatar
Mark Friedrichs committed
159
160
         delta    *= -1.0/rT2;
 
161
    for (int ii = 0; ii < 3; ii++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
162
163
164
       deltaR[P][ii]  = positions[particle2][ii] + deltaR[T][ii]*delta;
       deltaR[AP][ii] = positions[particle1][ii] - deltaR[P][ii];
       deltaR[CP][ii] = positions[particle3][ii] - deltaR[P][ii];
Jason Swails's avatar
More  
Jason Swails committed
165
    }
Mark Friedrichs's avatar
Mark Friedrichs committed
166
 
167
168
169
    double rAp2 = dotVector3(deltaR[AP],  deltaR[AP]);
    double rCp2 = dotVector3(deltaR[CP],  deltaR[CP]);
    if (rAp2 <= 0.0 && rCp2 <= 0.0) {
Mark Friedrichs's avatar
Mark Friedrichs committed
170
171
172
        return;
    }

173
    crossProductVector3(deltaR[CP], deltaR[AP], deltaR[M]);
Mark Friedrichs's avatar
Mark Friedrichs committed
174
 
175
176
177
    double rm = dotVector3(deltaR[M], deltaR[M]);
         rm   = sqrt(rm);
    if (rm < 0.000001) {
Mark Friedrichs's avatar
Mark Friedrichs committed
178
179
180
       rm = 0.000001;
    }
 
181
182
    double dot     = dotVector3(deltaR[AP], deltaR[CP]);
    double cosine  = dot/sqrt(rAp2*rCp2);
Mark Friedrichs's avatar
Mark Friedrichs committed
183
184
185
 
    double dEdR;
    double energyTerm;
186
    getPrefactorsGivenInPlaneAngleCosine(cosine, idealInPlaneAngle, quadraticK, cubicK,
peastman's avatar
peastman committed
187
                                         quarticK, penticK, sexticK, &dEdR,  &energyTerm);
Mark Friedrichs's avatar
Mark Friedrichs committed
188
189
190
191
 
    double termA   = -dEdR/(rAp2*rm);
    double termC   =  dEdR/(rCp2*rm);
 
192
193
    crossProductVector3(deltaR[AP], deltaR[M], deltaR[APxM]);
    crossProductVector3(deltaR[CP], deltaR[M], deltaR[CPxM]);
Mark Friedrichs's avatar
Mark Friedrichs committed
194
195
196
197
198
199
 
    // forces will be gathered here
 
    enum { dA, dB, dC, dD, LastDIndex };
    double forceTerm[LastDIndex][3];
 
200
    for (int ii = 0; ii < 3; ii++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
201
202
       forceTerm[dA][ii] = deltaR[APxM][ii]*termA;
       forceTerm[dC][ii] = deltaR[CPxM][ii]*termC;
203
       forceTerm[dB][ii] = -1.0*(forceTerm[dA][ii] + forceTerm[dC][ii]);
Mark Friedrichs's avatar
Mark Friedrichs committed
204
205
    }
 
206
    double pTrT2  = dotVector3(forceTerm[dB], deltaR[T]);
Mark Friedrichs's avatar
Mark Friedrichs committed
207
208
         pTrT2   /= rT2;
 
209
210
    crossProductVector3(deltaR[CD], forceTerm[dB], deltaR[CDxdB]);
    crossProductVector3(forceTerm[dB], deltaR[AD], deltaR[dBxAD]);
Mark Friedrichs's avatar
Mark Friedrichs committed
211
 
212
    if (fabs(pTrT2) > 1.0e-08) {
Mark Friedrichs's avatar
Mark Friedrichs committed
213
214
       double delta2 = delta*2.0;
 
215
       crossProductVector3(deltaR[BD], deltaR[CD], deltaR[BDxCD]);
216
       crossProductVector3(deltaR[T],  deltaR[CD], deltaR[TxCD]);
217
       crossProductVector3(deltaR[AD], deltaR[BD], deltaR[ADxBD]);
218
       crossProductVector3(deltaR[AD], deltaR[T],  deltaR[ADxT]);
219
       for (int ii = 0; ii < 3; ii++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
220
221
222
223
224
225
226
 
          double term           = deltaR[BDxCD][ii] + delta2*deltaR[TxCD][ii];
          forceTerm[dA][ii]  += delta*deltaR[CDxdB][ii] + term*pTrT2;
 
               term           = deltaR[ADxBD][ii] + delta2*deltaR[ADxT][ii];
          forceTerm[dC][ii]  += delta*deltaR[dBxAD][ii] + term*pTrT2;
 
227
          forceTerm[dD][ii]  = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
Mark Friedrichs's avatar
Mark Friedrichs committed
228
       }
Jason Swails's avatar
More  
Jason Swails committed
229
230
    }
    else {
231
       for (int ii = 0; ii < 3; ii++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
232
233
234
235
 
          forceTerm[dA][ii] += delta*deltaR[CDxdB][ii];
          forceTerm[dC][ii] += delta*deltaR[dBxAD][ii];
 
236
          forceTerm[dD][ii]  = -(forceTerm[dA][ii] + forceTerm[dB][ii] + forceTerm[dC][ii]);
Mark Friedrichs's avatar
Mark Friedrichs committed
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
       }
    }
 
    // accumulate forces and energy
 
    *energy                    += energyTerm;
 
    forces[particle1][0]       -= forceTerm[0][0];
    forces[particle1][1]       -= forceTerm[0][1];
    forces[particle1][2]       -= forceTerm[0][2];

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

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

    forces[particle4][0]       -= forceTerm[3][0];
    forces[particle4][1]       -= forceTerm[3][1];
    forces[particle4][2]       -= forceTerm[3][2];

}

262
static void computeAmoebaInPlaneAngleForces(Context& context, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce,
peastman's avatar
peastman committed
263
                                                     std::vector<Vec3>& expectedForces, double* expectedEnergy) {
Mark Friedrichs's avatar
Mark Friedrichs committed
264
265
266
267
268

    // get positions and zero forces

    State state                 = context.getState(State::Positions);
    std::vector<Vec3> positions = state.getPositions();
269
    expectedForces.resize(positions.size());
Mark Friedrichs's avatar
Mark Friedrichs committed
270
    
271
    for (unsigned int ii = 0; ii < expectedForces.size(); ii++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
272
273
274
275
276
277
        expectedForces[ii][0] = expectedForces[ii][1] = expectedForces[ii][2] = 0.0;
    }

    // calculates forces/energy

    *expectedEnergy = 0.0;
278
    for (int ii = 0; ii < AmoebaInPlaneAngleForce.getNumAngles(); ii++) {
peastman's avatar
peastman committed
279
        computeAmoebaInPlaneAngleForce(ii, positions, AmoebaInPlaneAngleForce, expectedForces, expectedEnergy);
Mark Friedrichs's avatar
Mark Friedrichs committed
280
281
282
    }
}

283
void compareWithExpectedForceAndEnergy(Context& context, AmoebaInPlaneAngleForce& AmoebaInPlaneAngleForce,
peastman's avatar
peastman committed
284
                                       double tolerance, const std::string& idString) {
Mark Friedrichs's avatar
Mark Friedrichs committed
285
286
287

    std::vector<Vec3> expectedForces;
    double expectedEnergy;
peastman's avatar
peastman committed
288
    computeAmoebaInPlaneAngleForces(context, AmoebaInPlaneAngleForce, expectedForces, &expectedEnergy);
Mark Friedrichs's avatar
Mark Friedrichs committed
289
290
291
292
   
    State state                      = context.getState(State::Forces | State::Energy);
    const std::vector<Vec3> forces   = state.getForces();

293
294
    for (unsigned int ii = 0; ii < forces.size(); ii++) {
        ASSERT_EQUAL_VEC(expectedForces[ii], forces[ii], tolerance);
Mark Friedrichs's avatar
Mark Friedrichs committed
295
    }
296
    ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), tolerance);
Mark Friedrichs's avatar
Mark Friedrichs committed
297
298
}

peastman's avatar
peastman committed
299
void testOneAngle() {
Mark Friedrichs's avatar
Mark Friedrichs committed
300
301
302

    System system;
    int numberOfParticles = 4;
303
    for (int ii = 0; ii < numberOfParticles; ii++) {
Mark Friedrichs's avatar
Mark Friedrichs committed
304
305
306
307
308
        system.addParticle(1.0);
    }

    LangevinIntegrator integrator(0.0, 0.1, 0.01);

309
    AmoebaInPlaneAngleForce* amoebaInPlaneAngleForce = new AmoebaInPlaneAngleForce();
Mark Friedrichs's avatar
Mark Friedrichs committed
310
311
312
313
314
315
316

    double angle      = 65.0;
    double quadraticK = 1.0;
    double cubicK     = 0.0e-01;
    double quarticK   = 0.0e-02;
    double penticK    = 0.0e-03;
    double sexticK    = 0.0e-04;
317
    amoebaInPlaneAngleForce->addAngle(0, 1, 2, 3, angle, quadraticK);
Mark Friedrichs's avatar
Mark Friedrichs committed
318

319
320
321
322
    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleCubic(cubicK);
    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleQuartic(quarticK);
    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAnglePentic(penticK);
    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleSextic(sexticK);
Mark Friedrichs's avatar
Mark Friedrichs committed
323

324
    system.addForce(amoebaInPlaneAngleForce);
325
326
    ASSERT(!amoebaInPlaneAngleForce->usesPeriodicBoundaryConditions());
    ASSERT(!system.usesPeriodicBoundaryConditions());
327
    Context context(system, integrator, Platform::getPlatformByName("Reference"));
Mark Friedrichs's avatar
Mark Friedrichs committed
328
329
330
331
332
333
334
335
336

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

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

    context.setPositions(positions);
peastman's avatar
peastman committed
337
    compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
338
339
340
341
342
343
344
    
    // Try changing the angle parameters and make sure it's still correct.
    
    amoebaInPlaneAngleForce->setAngleParameters(0, 0, 1, 2, 3, 1.1*angle, 1.4*quadraticK);
    bool exceptionThrown = false;
    try {
        // This should throw an exception.
peastman's avatar
peastman committed
345
        compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
346
347
348
349
350
351
    }
    catch (std::exception ex) {
        exceptionThrown = true;
    }
    ASSERT(exceptionThrown);
    amoebaInPlaneAngleForce->updateParametersInContext(context);
peastman's avatar
peastman committed
352
    compareWithExpectedForceAndEnergy(context, *amoebaInPlaneAngleForce, TOL, "testOneInPlaneAngle");
Mark Friedrichs's avatar
Mark Friedrichs committed
353
354
}

355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
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 = 4;
    for (int ii = 0; ii < numberOfParticles; ii++)
        system.addParticle(1.0);
    LangevinIntegrator integrator(0.0, 0.1, 0.01);
    AmoebaInPlaneAngleForce* amoebaInPlaneAngleForce = new AmoebaInPlaneAngleForce();
    double angle      = 65.0;
    double quadraticK = 1.0;
    double cubicK     = 0.0e-01;
    double quarticK   = 0.0e-02;
    double penticK    = 0.0e-03;
    double sexticK    = 0.0e-04;
    amoebaInPlaneAngleForce->addAngle(0, 1, 2, 3, angle, quadraticK);
    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleCubic(cubicK);
    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleQuartic(quarticK);
    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAnglePentic(penticK);
    amoebaInPlaneAngleForce->setAmoebaGlobalInPlaneAngleSextic(sexticK);
    amoebaInPlaneAngleForce->setUsesPeriodicBoundaryConditions(true);
    system.addForce(amoebaInPlaneAngleForce);
    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);
    positions[3] = Vec3(1, 1, 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);
}

399
int main(int numberOfArguments, char* argv[]) {
Mark Friedrichs's avatar
Mark Friedrichs committed
400
401

    try {
402
        std::cout << "TestReferenceAmoebaInPlaneAngleForce running test..." << std::endl;
403
        registerAmoebaReferenceKernelFactories();
peastman's avatar
peastman committed
404
        testOneAngle();
405
        testPeriodic();
Mark Friedrichs's avatar
Mark Friedrichs committed
406
407
408
409
410
411
    }
    catch(const std::exception& e) {
        std::cout << "exception: " << e.what() << std::endl;
        std::cout << "FAIL - ERROR.  Test failed." << std::endl;
        return 1;
    }
412
413
    //std::cout << "PASS - Test succeeded." << std::endl;
    std::cout << "Done" << std::endl;
Mark Friedrichs's avatar
Mark Friedrichs committed
414
415
    return 0;
}