"vscode:/vscode.git/clone" did not exist on "1f6f23ceb054d8f40039d41cecb7ad8677b0d694"
Commit 68ef42c0 authored by Rossen Apostolov's avatar Rossen Apostolov
Browse files

New tests for Ewald and PME on both Reference and Cuda platforms.

parent 25af59f8
/* -------------------------------------------------------------------------- *
* 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. *
* -------------------------------------------------------------------------- */
/**
* This tests the Ewald summation method cuda implementation of NonbondedForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "ReferencePlatform.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/internal/ContextImpl.h"
#include "kernels/gputypes.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include "../src/sfmt/SFMT.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testEwaldPME() {
// Use amorphous NaCl system for the tests
const int numParticles = 216;
const double cutoff = 0.8;
const double boxSize = 1.86206;
const double tol = 0.001;
CudaPlatform cuda;
ReferencePlatform reference;
System system;
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(tol);
for (int i = 0; i < numParticles/2; i++)
system.addParticle(22.99);
for (int i = 0; i < numParticles/2; i++)
system.addParticle(35.45);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(1.0, 1.0,0.0);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(-1.0, 1.0,0.0);
system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
system.addForce(nonbonded);
vector<Vec3> positions(numParticles);
#include "nacl_amorph.dat"
// (1) Check whether the Reference and Cuda platforms agree when using Ewald Method
Context cudaContext(system, integrator, cuda);
Context referenceContext(system, integrator, reference);
cudaContext.setPositions(positions);
referenceContext.setPositions(positions);
State cudaState = cudaContext.getState(State::Forces | State::Energy);
State referenceState = referenceContext.getState(State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(cudaState.getForces()[i], referenceState.getForces()[i], tol);
}
ASSERT_EQUAL_TOL(cudaState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
// (2) Check whether Ewald method in Cuda is self-consistent
double norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = cudaState.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-3;
double step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = cudaState.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
Context cudaContext2(system, integrator, cuda);
cudaContext2.setPositions(positions);
State cudaState2 = cudaContext2.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (cudaState2.getPotentialEnergy()-cudaState.getPotentialEnergy())/delta, tol)
// (3) Check whether the Reference and Cuda platforms agree when using PME
nonbonded->setNonbondedMethod(NonbondedForce::PME);
cudaContext.reinitialize();
referenceContext.reinitialize();
cudaContext.setPositions(positions);
referenceContext.setPositions(positions);
cudaState = cudaContext.getState(State::Forces | State::Energy);
referenceState = referenceContext.getState(State::Forces | State::Energy);
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(cudaState.getForces()[i], referenceState.getForces()[i], tol);
}
ASSERT_EQUAL_TOL(cudaState.getPotentialEnergy(), referenceState.getPotentialEnergy(), tol);
// (4) Check whether PME method in Cuda is self-consistent
norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = cudaState.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = cudaState.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
Context cudaContext3(system, integrator, cuda);
cudaContext3.setPositions(positions);
State cudaState3 = cudaContext3.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (cudaState3.getPotentialEnergy()-cudaState.getPotentialEnergy())/delta, tol)
}
void testEwald2Ions() {
CudaPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(-1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(TOL);
system.setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(3.048000,2.764000,3.156000);
positions[1] = Vec3(2.809000,2.888000,2.571000);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
ASSERT_EQUAL_VEC(Vec3( 123.711, -64.1877, 302.716), forces[1], 10*TOL);
ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 0.01/*10*TOL*/);
}
int main() {
try {
testEwaldPME();
// testEwald2Ions();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
......@@ -359,34 +359,6 @@ void testPeriodic() {
ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
}
void testEwald() {
CudaPlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(-1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(TOL);
system.setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(3.048000,2.764000,3.156000);
positions[1] = Vec3(2.809000,2.888000,2.571000);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
ASSERT_EQUAL_VEC(Vec3( 123.711, -64.1877, 302.716), forces[1], 10*TOL);
ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 0.01/*10*TOL*/);
}
void testLargeSystem() {
const int numMolecules = 600;
......@@ -636,7 +608,6 @@ int main() {
testCutoff();
testCutoff14();
testPeriodic();
testEwald();
testLargeSystem();
testBlockInteractions(false);
testBlockInteractions(true);
......
positions[0] = Vec3(0.230000,0.628000,0.113000);
positions[1] = Vec3(0.225000,0.275000,-0.866000);
positions[2] = Vec3(0.019000,0.368000,0.647000);
positions[3] = Vec3(0.569000,-0.587000,-0.697000);
positions[4] = Vec3(-0.307000,-0.351000,0.703000);
positions[5] = Vec3(-0.119000,0.618000,0.856000);
positions[6] = Vec3(-0.727000,0.703000,0.717000);
positions[7] = Vec3(-0.107000,0.607000,0.231000);
positions[8] = Vec3(0.768000,-0.718000,-0.839000);
positions[9] = Vec3(0.850000,0.798000,-0.039000);
positions[10] = Vec3(0.685000,-0.850000,0.665000);
positions[11] = Vec3(0.686000,-0.701000,-0.059000);
positions[12] = Vec3(0.335000,-0.427000,-0.801000);
positions[13] = Vec3(-0.402000,-0.357000,-0.523000);
positions[14] = Vec3(0.438000,0.392000,-0.363000);
positions[15] = Vec3(-0.259000,0.447000,0.737000);
positions[16] = Vec3(0.231000,-0.149000,0.483000);
positions[17] = Vec3(-0.735000,-0.521000,-0.172000);
positions[18] = Vec3(0.230000,-0.428000,0.538000);
positions[19] = Vec3(0.240000,-0.771000,0.886000);
positions[20] = Vec3(0.620000,-0.076000,-0.423000);
positions[21] = Vec3(0.606000,-0.898000,0.123000);
positions[22] = Vec3(-0.268000,0.114000,-0.382000);
positions[23] = Vec3(0.122000,0.643000,0.563000);
positions[24] = Vec3(-0.020000,-0.095000,0.359000);
positions[25] = Vec3(0.027000,-0.266000,0.117000);
positions[26] = Vec3(-0.173000,0.922000,0.612000);
positions[27] = Vec3(-0.221000,-0.754000,0.432000);
positions[28] = Vec3(0.113000,0.737000,-0.265000);
positions[29] = Vec3(0.613000,-0.497000,0.726000);
positions[30] = Vec3(-0.569000,-0.634000,-0.439000);
positions[31] = Vec3(0.809000,0.004000,0.502000);
positions[32] = Vec3(0.197000,-0.886000,-0.598000);
positions[33] = Vec3(-0.337000,-0.863000,0.190000);
positions[34] = Vec3(-0.675000,-0.070000,-0.246000);
positions[35] = Vec3(0.317000,0.251000,-0.061000);
positions[36] = Vec3(-0.396000,-0.445000,-0.909000);
positions[37] = Vec3(-0.195000,-0.148000,0.572000);
positions[38] = Vec3(0.598000,0.729000,0.270000);
positions[39] = Vec3(-0.581000,0.345000,-0.918000);
positions[40] = Vec3(-0.286000,-0.200000,0.307000);
positions[41] = Vec3(0.807000,0.605000,-0.397000);
positions[42] = Vec3(-0.468000,0.469000,-0.188000);
positions[43] = Vec3(-0.889000,0.890000,-0.290000);
positions[44] = Vec3(-0.871000,0.410000,-0.620000);
positions[45] = Vec3(-0.821000,0.701000,0.429000);
positions[46] = Vec3(0.076000,0.811000,0.789000);
positions[47] = Vec3(0.130000,-0.041000,-0.291000);
positions[48] = Vec3(0.865000,0.348000,0.195000);
positions[49] = Vec3(-0.143000,0.585000,-0.031000);
positions[50] = Vec3(-0.500000,-0.718000,0.545000);
positions[51] = Vec3(0.550000,0.196000,0.885000);
positions[52] = Vec3(-0.854000,-0.406000,0.477000);
positions[53] = Vec3(0.351000,-0.061000,0.853000);
positions[54] = Vec3(-0.067000,-0.796000,0.873000);
positions[55] = Vec3(-0.635000,-0.312000,-0.356000);
positions[56] = Vec3(0.321000,-0.919000,0.242000);
positions[57] = Vec3(-0.404000,0.735000,0.728000);
positions[58] = Vec3(0.461000,-0.596000,-0.135000);
positions[59] = Vec3(-0.751000,-0.086000,0.237000);
positions[60] = Vec3(0.202000,0.285000,-0.364000);
positions[61] = Vec3(-0.230000,-0.485000,0.081000);
positions[62] = Vec3(0.464000,-0.119000,0.323000);
positions[63] = Vec3(-0.462000,0.107000,0.426000);
positions[64] = Vec3(0.249000,-0.077000,-0.621000);
positions[65] = Vec3(-0.922000,-0.164000,0.904000);
positions[66] = Vec3(0.382000,0.700000,0.480000);
positions[67] = Vec3(-0.315000,0.222000,-0.133000);
positions[68] = Vec3(0.614000,0.122000,0.117000);
positions[69] = Vec3(0.781000,0.264000,-0.113000);
positions[70] = Vec3(0.888000,-0.348000,-0.667000);
positions[71] = Vec3(-0.511000,0.590000,-0.429000);
positions[72] = Vec3(0.803000,-0.460000,0.924000);
positions[73] = Vec3(0.922000,0.503000,0.899000);
positions[74] = Vec3(0.539000,0.064000,0.512000);
positions[75] = Vec3(-0.428000,-0.674000,0.041000);
positions[76] = Vec3(0.297000,0.035000,0.171000);
positions[77] = Vec3(-0.927000,0.236000,0.480000);
positions[78] = Vec3(-0.786000,0.683000,-0.398000);
positions[79] = Vec3(-0.635000,-0.292000,0.793000);
positions[80] = Vec3(0.459000,-0.710000,0.741000);
positions[81] = Vec3(-0.591000,-0.065000,0.591000);
positions[82] = Vec3(-0.830000,0.549000,0.016000);
positions[83] = Vec3(0.078000,0.556000,-0.476000);
positions[84] = Vec3(0.561000,0.222000,-0.715000);
positions[85] = Vec3(0.866000,0.454000,0.642000);
positions[86] = Vec3(-0.845000,0.039000,0.753000);
positions[87] = Vec3(-0.433000,-0.689000,0.867000);
positions[88] = Vec3(-0.396000,0.590000,-0.870000);
positions[89] = Vec3(-0.005000,0.833000,0.377000);
positions[90] = Vec3(0.488000,-0.477000,0.174000);
positions[91] = Vec3(-0.198000,-0.582000,0.657000);
positions[92] = Vec3(-0.472000,0.575000,0.078000);
positions[93] = Vec3(0.527000,0.256000,0.328000);
positions[94] = Vec3(-0.108000,-0.639000,-0.274000);
positions[95] = Vec3(-0.798000,-0.515000,-0.522000);
positions[96] = Vec3(-0.270000,-0.233000,-0.237000);
positions[97] = Vec3(-0.751000,-0.667000,-0.762000);
positions[98] = Vec3(-0.224000,-0.763000,-0.783000);
positions[99] = Vec3(0.915000,0.089000,-0.460000);
positions[100] = Vec3(-0.882000,-0.746000,-0.143000);
positions[101] = Vec3(0.705000,-0.812000,0.368000);
positions[102] = Vec3(0.410000,0.813000,-0.611000);
positions[103] = Vec3(-0.588000,0.386000,-0.600000);
positions[104] = Vec3(0.064000,-0.298000,-0.531000);
positions[105] = Vec3(0.367000,-0.762000,0.501000);
positions[106] = Vec3(0.566000,0.537000,0.865000);
positions[107] = Vec3(-0.610000,-0.514000,0.388000);
positions[108] = Vec3(-0.590000,-0.417000,-0.720000);
positions[109] = Vec3(-0.280000,0.639000,0.472000);
positions[110] = Vec3(0.354000,-0.352000,-0.533000);
positions[111] = Vec3(0.402000,0.751000,-0.264000);
positions[112] = Vec3(-0.275000,0.779000,-0.192000);
positions[113] = Vec3(-0.849000,0.105000,-0.092000);
positions[114] = Vec3(0.504000,0.050000,-0.122000);
positions[115] = Vec3(0.573000,0.870000,-0.833000);
positions[116] = Vec3(-0.502000,0.862000,-0.817000);
positions[117] = Vec3(-0.653000,0.525000,0.275000);
positions[118] = Vec3(0.307000,0.213000,-0.631000);
positions[119] = Vec3(0.037000,-0.552000,-0.580000);
positions[120] = Vec3(0.732000,0.634000,-0.798000);
positions[121] = Vec3(-0.134000,-0.927000,-0.008000);
positions[122] = Vec3(0.307000,0.063000,0.618000);
positions[123] = Vec3(-0.240000,0.367000,0.374000);
positions[124] = Vec3(-0.839000,0.766000,-0.896000);
positions[125] = Vec3(-0.882000,-0.289000,-0.162000);
positions[126] = Vec3(-0.003000,-0.344000,-0.257000);
positions[127] = Vec3(0.350000,0.898000,-0.058000);
positions[128] = Vec3(-0.322000,0.274000,0.125000);
positions[129] = Vec3(-0.559000,0.838000,0.042000);
positions[130] = Vec3(-0.794000,-0.529000,0.849000);
positions[131] = Vec3(0.319000,0.810000,-0.913000);
positions[132] = Vec3(0.339000,0.509000,-0.856000);
positions[133] = Vec3(0.511000,0.415000,-0.054000);
positions[134] = Vec3(-0.724000,0.380000,-0.184000);
positions[135] = Vec3(-0.702000,0.207000,-0.385000);
positions[136] = Vec3(0.008000,-0.536000,0.200000);
positions[137] = Vec3(0.088000,-0.061000,0.927000);
positions[138] = Vec3(0.504000,-0.294000,0.910000);
positions[139] = Vec3(-0.860000,0.796000,-0.624000);
positions[140] = Vec3(0.040000,0.544000,-0.748000);
positions[141] = Vec3(0.189000,0.520000,-0.140000);
positions[142] = Vec3(-0.493000,-0.912000,-0.202000);
positions[143] = Vec3(0.815000,0.572000,0.325000);
positions[144] = Vec3(-0.205000,0.604000,-0.656000);
positions[145] = Vec3(0.252000,-0.298000,-0.118000);
positions[146] = Vec3(0.671000,0.464000,-0.593000);
positions[147] = Vec3(0.930000,-0.184000,-0.397000);
positions[148] = Vec3(0.473000,0.500000,0.191000);
positions[149] = Vec3(0.159000,-0.725000,-0.396000);
positions[150] = Vec3(-0.515000,-0.803000,-0.628000);
positions[151] = Vec3(-0.560000,0.855000,0.309000);
positions[152] = Vec3(-0.103000,-0.115000,-0.708000);
positions[153] = Vec3(-0.610000,-0.131000,-0.734000);
positions[154] = Vec3(0.083000,-0.604000,-0.840000);
positions[155] = Vec3(0.688000,-0.200000,-0.146000);
positions[156] = Vec3(0.903000,0.086000,0.133000);
positions[157] = Vec3(-0.136000,0.135000,0.523000);
positions[158] = Vec3(-0.474000,-0.289000,0.477000);
positions[159] = Vec3(0.130000,-0.068000,-0.011000);
positions[160] = Vec3(-0.582000,0.927000,0.672000);
positions[161] = Vec3(0.830000,-0.589000,-0.440000);
positions[162] = Vec3(0.672000,-0.246000,0.154000);
positions[163] = Vec3(-0.212000,-0.142000,-0.468000);
positions[164] = Vec3(-0.021000,0.175000,-0.899000);
positions[165] = Vec3(0.263000,0.326000,0.720000);
positions[166] = Vec3(-0.668000,-0.250000,0.031000);
positions[167] = Vec3(0.822000,-0.860000,-0.490000);
positions[168] = Vec3(0.916000,0.910000,0.291000);
positions[169] = Vec3(-0.358000,-0.255000,0.044000);
positions[170] = Vec3(0.372000,-0.574000,-0.372000);
positions[171] = Vec3(-0.248000,-0.570000,-0.573000);
positions[172] = Vec3(-0.823000,-0.764000,0.696000);
positions[173] = Vec3(-0.848000,0.236000,-0.891000);
positions[174] = Vec3(0.590000,-0.375000,0.491000);
positions[175] = Vec3(-0.153000,0.385000,-0.481000);
positions[176] = Vec3(0.255000,-0.514000,0.290000);
positions[177] = Vec3(0.105000,-0.849000,-0.136000);
positions[178] = Vec3(0.672000,0.203000,-0.373000);
positions[179] = Vec3(0.075000,0.345000,0.033000);
positions[180] = Vec3(-0.422000,0.856000,-0.464000);
positions[181] = Vec3(0.072000,0.166000,0.318000);
positions[182] = Vec3(-0.679000,-0.527000,0.119000);
positions[183] = Vec3(0.613000,0.842000,-0.431000);
positions[184] = Vec3(-0.369000,-0.095000,-0.903000);
positions[185] = Vec3(0.716000,0.565000,-0.154000);
positions[186] = Vec3(-0.412000,-0.642000,-0.229000);
positions[187] = Vec3(0.390000,-0.121000,-0.302000);
positions[188] = Vec3(-0.188000,0.883000,-0.608000);
positions[189] = Vec3(-0.637000,0.325000,0.449000);
positions[190] = Vec3(0.594000,0.745000,0.652000);
positions[191] = Vec3(-0.085000,0.342000,-0.220000);
positions[192] = Vec3(-0.132000,-0.928000,-0.345000);
positions[193] = Vec3(0.859000,-0.488000,0.016000);
positions[194] = Vec3(0.661000,-0.072000,-0.909000);
positions[195] = Vec3(-0.454000,-0.011000,-0.142000);
positions[196] = Vec3(0.859000,-0.906000,0.861000);
positions[197] = Vec3(-0.779000,-0.878000,0.087000);
positions[198] = Vec3(-0.001000,-0.293000,0.851000);
positions[199] = Vec3(0.221000,-0.548000,-0.018000);
positions[200] = Vec3(0.079000,-0.622000,0.653000);
positions[201] = Vec3(0.672000,-0.471000,-0.238000);
positions[202] = Vec3(-0.038000,0.192000,-0.635000);
positions[203] = Vec3(0.428000,0.424000,0.520000);
positions[204] = Vec3(-0.157000,-0.375000,-0.758000);
positions[205] = Vec3(0.317000,0.547000,-0.582000);
positions[206] = Vec3(0.812000,-0.276000,0.687000);
positions[207] = Vec3(-0.438000,0.214000,-0.750000);
positions[208] = Vec3(-0.861000,0.034000,-0.708000);
positions[209] = Vec3(0.770000,-0.532000,0.301000);
positions[210] = Vec3(0.618000,-0.295000,-0.578000);
positions[211] = Vec3(-0.510000,0.052000,0.168000);
positions[212] = Vec3(-0.562000,0.453000,0.691000);
positions[213] = Vec3(-0.269000,0.221000,0.882000);
positions[214] = Vec3(0.039000,-0.785000,0.300000);
positions[215] = Vec3(0.875000,-0.216000,0.337000);
......@@ -49,42 +49,47 @@ using namespace std;
const double TOL = 1e-5;
void testLargeSystem() {
void testEwaldExact() {
// Use a NaCl crystal to compare the calculated and Madelung energies
const int numParticles = 1000;
const double cutoff = 1.0;
const double boxSize = 2.82;
ReferencePlatform platform;
System system;
for (int i = 0; i < 500; i++)
for (int i = 0; i < numParticles/2; i++)
system.addParticle(22.99);
for (int i = 0; i < 500; i++)
for (int i = 0; i < numParticles/2; i++)
system.addParticle(35.45);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < 500; i++)
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(1.0, 1.0,0.0);
// nonbonded->addParticle(1.0, 0.33284,0.0115897);
for (int i = 0; i < 500; i++)
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(-1.0, 1.0,0.0);
// nonbonded->addParticle(-1.0, 0.440104,0.4184);
nonbonded->setNonbondedMethod(NonbondedForce::PME);
const double cutoff = 1.0;
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
nonbonded->setCutoffDistance(cutoff);
system.setPeriodicBoxVectors(Vec3(2.82, 0, 0), Vec3(0, 2.82, 0), Vec3(0, 0, 2.82));
system.setPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
nonbonded->setEwaldErrorTolerance(TOL);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(1000);
vector<Vec3> positions(numParticles);
#include "nacl_crystal.dat"
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
// cout << "force 0: " << forces[0] << endl;
// cout << "force 1: " << forces[1] << endl;
cout << "PotentialEnergy: " << state.getPotentialEnergy() << endl;
ASSERT_EQUAL_TOL(-430355.0, state.getPotentialEnergy(), 100*TOL);
// ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
// ASSERT_EQUAL_VEC(Vec3(123.711, -64.1877, 302.716), forces[1], 10*TOL);
// Compare against Gromacs result for now
ASSERT_EQUAL_TOL(-430494.0, state.getPotentialEnergy(), 10*TOL);
}
void testEwald() {
void testEwald2Ions() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
......@@ -93,9 +98,6 @@ void testEwald() {
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(-1.0, 1, 0);
// Sodium Chloride
// nonbonded->addParticle(1.0, 0.33284,0.0115897);
// nonbonded->addParticle(-1.0, 0.440104,0.4184);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
......@@ -109,53 +111,184 @@ void testEwald() {
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
// cout << "force 0: " << forces[0] << endl;
// cout << "force 1: " << forces[1] << endl;
// cout << "PotentialEnergy: " << state.getPotentialEnergy() << endl;
ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
ASSERT_EQUAL_VEC(Vec3(123.711, -64.1877, 302.716), forces[1], 10*TOL);
ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 10*TOL);
}
void testPME() {
void testWaterSystem() {
ReferencePlatform platform;
System system;
for (int i = 0 ; i < 42 ; i++)
static int numParticles;
numParticles = 648;
for (int i = 0 ; i < numParticles ; i++)
{
system.addParticle(1.0);
}
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0 ; i < 14 ; i++)
for (int i = 0 ; i < numParticles/3 ; i++)
{
nonbonded->addParticle(-0.82, 1, 0);
nonbonded->addParticle(0.41, 1, 0);
nonbonded->addParticle(0.41, 1, 0);
}
nonbonded->setNonbondedMethod(NonbondedForce::PME);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 0.8;
nonbonded->setCutoffDistance(cutoff);
system.setPeriodicBoxVectors(Vec3(1.86206, 0, 0), Vec3(0, 1.86206, 0), Vec3(0, 0, 1.86206));
nonbonded->setEwaldErrorTolerance(TOL);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(42);
vector<Vec3> positions(numParticles);
#include "water.dat"
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
State state1 = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state1.getForces();
// for (int i = 0 ; i < 42 ; i++)
// cout << "f [" << i << " : ]" << forces[i] << endl;
// cout << "PotentialEnergy: " << state.getPotentialEnergy() << endl;
// ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
// ASSERT_EQUAL_VEC(Vec3(123.711, -64.1877, 302.716), forces[1], 10*TOL);
// cout << "PotentialEnergy: " << state1.getPotentialEnergy() << endl;
// Take a small step in the direction of the energy gradient.
double norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = state1.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-3;
double step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = state1.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
context.setPositions(positions);
// See whether the potential energy changed by the expected amount.
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
State state2 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state1.getPotentialEnergy())/delta, 0.01)
}
void testEwaldPME() {
// Use amorphoush NaCl system
ReferencePlatform platform;
const int numParticles = 216;
System system;
for (int i = 0; i < numParticles/2; i++)
system.addParticle(22.99);
for (int i = 0; i < numParticles/2; i++)
system.addParticle(35.45);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(1.0, 1.0,0.0);
for (int i = 0; i < numParticles/2; i++)
nonbonded->addParticle(-1.0, 1.0,0.0);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 0.8;
nonbonded->setCutoffDistance(cutoff);
system.setPeriodicBoxVectors(Vec3(1.86206, 0, 0), Vec3(0, 1.86206, 0), Vec3(0, 0, 1.86206));
nonbonded->setEwaldErrorTolerance(TOL);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
#include "nacl_amorph.dat"
context.setPositions(positions);
State state1 = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces1 = state1.getForces();
// (1) CHECK EXACT VALUE OF EWALD ENERGY
ASSERT_EQUAL_TOL(-26651.9, state1.getPotentialEnergy(), TOL);
// (2) CHECK WHETHER THE EWALD FORCES ARE THE SAME AS THE GROMACS OUTPUT
// Even at tolerance 0.1 the test doesn't pass
// #include "nacl_amorph_GromacsForcesEwald.dat"
// (3) CHECK SELF-CONSISTENCY
// Take a small step in the direction of the energy gradient.
double norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = state1.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-3;
double step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = state1.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
context.setPositions(positions);
// See whether the potential energy changed by the expected amount.
State state2 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state1.getPotentialEnergy())/delta, 0.01)
// (4) CHECK EXACT VALUE OF PME ENERGY
nonbonded->setNonbondedMethod(NonbondedForce::PME);
context.reinitialize();
#include "nacl_amorph.dat"
context.setPositions(positions);
State state3 = context.getState(State::Forces | State::Energy);
ASSERT_EQUAL_TOL(-26651.9, state3.getPotentialEnergy(), 10*TOL);
// (5) CHECK WHETHER PME FORCES ARE THE SAME AS THE GROMACS OUTPUT
// Even at tolerance 0.1 the test doesn't pass
// #include "nacl_amorph_GromacsForcesEwald.dat"
// (6) CHECK PME FOR SELF-CONSISTENCY
// Take a small step in the direction of the energy gradient.
norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = state3.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = state3.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
context.setPositions(positions);
// See whether the potential energy changed by the expected amount.
State state4 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state4.getPotentialEnergy()-state3.getPotentialEnergy())/delta, 0.01)
}
int main() {
try {
testLargeSystem();
testEwald();
testPME();
testEwaldExact();
// testEwald2Ions();
// testWaterSystem();
testEwaldPME();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
......@@ -355,34 +355,6 @@ void testPeriodic() {
ASSERT_EQUAL_TOL(2*138.935485*(1.0)*(1.0+krf*1.0-crf), state.getPotentialEnergy(), TOL);
}
void testEwald() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(-1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
nonbonded->setEwaldErrorTolerance(TOL);
system.setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(3.048000,2.764000,3.156000);
positions[1] = Vec3(2.809000,2.888000,2.571000);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
ASSERT_EQUAL_VEC(Vec3( 123.711, -64.1877, 302.716), forces[1], 10*TOL);
ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 0.01/*10*TOL*/);
}
int main() {
try {
testCoulomb();
......@@ -391,7 +363,6 @@ int main() {
testCutoff();
testCutoff14();
testPeriodic();
testEwald();
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
positions[0] = Vec3(0.230000,0.628000,0.113000);
positions[1] = Vec3(0.225000,0.275000,-0.866000);
positions[2] = Vec3(0.019000,0.368000,0.647000);
positions[3] = Vec3(0.569000,-0.587000,-0.697000);
positions[4] = Vec3(-0.307000,-0.351000,0.703000);
positions[5] = Vec3(-0.119000,0.618000,0.856000);
positions[6] = Vec3(-0.727000,0.703000,0.717000);
positions[7] = Vec3(-0.107000,0.607000,0.231000);
positions[8] = Vec3(0.768000,-0.718000,-0.839000);
positions[9] = Vec3(0.850000,0.798000,-0.039000);
positions[10] = Vec3(0.685000,-0.850000,0.665000);
positions[11] = Vec3(0.686000,-0.701000,-0.059000);
positions[12] = Vec3(0.335000,-0.427000,-0.801000);
positions[13] = Vec3(-0.402000,-0.357000,-0.523000);
positions[14] = Vec3(0.438000,0.392000,-0.363000);
positions[15] = Vec3(-0.259000,0.447000,0.737000);
positions[16] = Vec3(0.231000,-0.149000,0.483000);
positions[17] = Vec3(-0.735000,-0.521000,-0.172000);
positions[18] = Vec3(0.230000,-0.428000,0.538000);
positions[19] = Vec3(0.240000,-0.771000,0.886000);
positions[20] = Vec3(0.620000,-0.076000,-0.423000);
positions[21] = Vec3(0.606000,-0.898000,0.123000);
positions[22] = Vec3(-0.268000,0.114000,-0.382000);
positions[23] = Vec3(0.122000,0.643000,0.563000);
positions[24] = Vec3(-0.020000,-0.095000,0.359000);
positions[25] = Vec3(0.027000,-0.266000,0.117000);
positions[26] = Vec3(-0.173000,0.922000,0.612000);
positions[27] = Vec3(-0.221000,-0.754000,0.432000);
positions[28] = Vec3(0.113000,0.737000,-0.265000);
positions[29] = Vec3(0.613000,-0.497000,0.726000);
positions[30] = Vec3(-0.569000,-0.634000,-0.439000);
positions[31] = Vec3(0.809000,0.004000,0.502000);
positions[32] = Vec3(0.197000,-0.886000,-0.598000);
positions[33] = Vec3(-0.337000,-0.863000,0.190000);
positions[34] = Vec3(-0.675000,-0.070000,-0.246000);
positions[35] = Vec3(0.317000,0.251000,-0.061000);
positions[36] = Vec3(-0.396000,-0.445000,-0.909000);
positions[37] = Vec3(-0.195000,-0.148000,0.572000);
positions[38] = Vec3(0.598000,0.729000,0.270000);
positions[39] = Vec3(-0.581000,0.345000,-0.918000);
positions[40] = Vec3(-0.286000,-0.200000,0.307000);
positions[41] = Vec3(0.807000,0.605000,-0.397000);
positions[42] = Vec3(-0.468000,0.469000,-0.188000);
positions[43] = Vec3(-0.889000,0.890000,-0.290000);
positions[44] = Vec3(-0.871000,0.410000,-0.620000);
positions[45] = Vec3(-0.821000,0.701000,0.429000);
positions[46] = Vec3(0.076000,0.811000,0.789000);
positions[47] = Vec3(0.130000,-0.041000,-0.291000);
positions[48] = Vec3(0.865000,0.348000,0.195000);
positions[49] = Vec3(-0.143000,0.585000,-0.031000);
positions[50] = Vec3(-0.500000,-0.718000,0.545000);
positions[51] = Vec3(0.550000,0.196000,0.885000);
positions[52] = Vec3(-0.854000,-0.406000,0.477000);
positions[53] = Vec3(0.351000,-0.061000,0.853000);
positions[54] = Vec3(-0.067000,-0.796000,0.873000);
positions[55] = Vec3(-0.635000,-0.312000,-0.356000);
positions[56] = Vec3(0.321000,-0.919000,0.242000);
positions[57] = Vec3(-0.404000,0.735000,0.728000);
positions[58] = Vec3(0.461000,-0.596000,-0.135000);
positions[59] = Vec3(-0.751000,-0.086000,0.237000);
positions[60] = Vec3(0.202000,0.285000,-0.364000);
positions[61] = Vec3(-0.230000,-0.485000,0.081000);
positions[62] = Vec3(0.464000,-0.119000,0.323000);
positions[63] = Vec3(-0.462000,0.107000,0.426000);
positions[64] = Vec3(0.249000,-0.077000,-0.621000);
positions[65] = Vec3(-0.922000,-0.164000,0.904000);
positions[66] = Vec3(0.382000,0.700000,0.480000);
positions[67] = Vec3(-0.315000,0.222000,-0.133000);
positions[68] = Vec3(0.614000,0.122000,0.117000);
positions[69] = Vec3(0.781000,0.264000,-0.113000);
positions[70] = Vec3(0.888000,-0.348000,-0.667000);
positions[71] = Vec3(-0.511000,0.590000,-0.429000);
positions[72] = Vec3(0.803000,-0.460000,0.924000);
positions[73] = Vec3(0.922000,0.503000,0.899000);
positions[74] = Vec3(0.539000,0.064000,0.512000);
positions[75] = Vec3(-0.428000,-0.674000,0.041000);
positions[76] = Vec3(0.297000,0.035000,0.171000);
positions[77] = Vec3(-0.927000,0.236000,0.480000);
positions[78] = Vec3(-0.786000,0.683000,-0.398000);
positions[79] = Vec3(-0.635000,-0.292000,0.793000);
positions[80] = Vec3(0.459000,-0.710000,0.741000);
positions[81] = Vec3(-0.591000,-0.065000,0.591000);
positions[82] = Vec3(-0.830000,0.549000,0.016000);
positions[83] = Vec3(0.078000,0.556000,-0.476000);
positions[84] = Vec3(0.561000,0.222000,-0.715000);
positions[85] = Vec3(0.866000,0.454000,0.642000);
positions[86] = Vec3(-0.845000,0.039000,0.753000);
positions[87] = Vec3(-0.433000,-0.689000,0.867000);
positions[88] = Vec3(-0.396000,0.590000,-0.870000);
positions[89] = Vec3(-0.005000,0.833000,0.377000);
positions[90] = Vec3(0.488000,-0.477000,0.174000);
positions[91] = Vec3(-0.198000,-0.582000,0.657000);
positions[92] = Vec3(-0.472000,0.575000,0.078000);
positions[93] = Vec3(0.527000,0.256000,0.328000);
positions[94] = Vec3(-0.108000,-0.639000,-0.274000);
positions[95] = Vec3(-0.798000,-0.515000,-0.522000);
positions[96] = Vec3(-0.270000,-0.233000,-0.237000);
positions[97] = Vec3(-0.751000,-0.667000,-0.762000);
positions[98] = Vec3(-0.224000,-0.763000,-0.783000);
positions[99] = Vec3(0.915000,0.089000,-0.460000);
positions[100] = Vec3(-0.882000,-0.746000,-0.143000);
positions[101] = Vec3(0.705000,-0.812000,0.368000);
positions[102] = Vec3(0.410000,0.813000,-0.611000);
positions[103] = Vec3(-0.588000,0.386000,-0.600000);
positions[104] = Vec3(0.064000,-0.298000,-0.531000);
positions[105] = Vec3(0.367000,-0.762000,0.501000);
positions[106] = Vec3(0.566000,0.537000,0.865000);
positions[107] = Vec3(-0.610000,-0.514000,0.388000);
positions[108] = Vec3(-0.590000,-0.417000,-0.720000);
positions[109] = Vec3(-0.280000,0.639000,0.472000);
positions[110] = Vec3(0.354000,-0.352000,-0.533000);
positions[111] = Vec3(0.402000,0.751000,-0.264000);
positions[112] = Vec3(-0.275000,0.779000,-0.192000);
positions[113] = Vec3(-0.849000,0.105000,-0.092000);
positions[114] = Vec3(0.504000,0.050000,-0.122000);
positions[115] = Vec3(0.573000,0.870000,-0.833000);
positions[116] = Vec3(-0.502000,0.862000,-0.817000);
positions[117] = Vec3(-0.653000,0.525000,0.275000);
positions[118] = Vec3(0.307000,0.213000,-0.631000);
positions[119] = Vec3(0.037000,-0.552000,-0.580000);
positions[120] = Vec3(0.732000,0.634000,-0.798000);
positions[121] = Vec3(-0.134000,-0.927000,-0.008000);
positions[122] = Vec3(0.307000,0.063000,0.618000);
positions[123] = Vec3(-0.240000,0.367000,0.374000);
positions[124] = Vec3(-0.839000,0.766000,-0.896000);
positions[125] = Vec3(-0.882000,-0.289000,-0.162000);
positions[126] = Vec3(-0.003000,-0.344000,-0.257000);
positions[127] = Vec3(0.350000,0.898000,-0.058000);
positions[128] = Vec3(-0.322000,0.274000,0.125000);
positions[129] = Vec3(-0.559000,0.838000,0.042000);
positions[130] = Vec3(-0.794000,-0.529000,0.849000);
positions[131] = Vec3(0.319000,0.810000,-0.913000);
positions[132] = Vec3(0.339000,0.509000,-0.856000);
positions[133] = Vec3(0.511000,0.415000,-0.054000);
positions[134] = Vec3(-0.724000,0.380000,-0.184000);
positions[135] = Vec3(-0.702000,0.207000,-0.385000);
positions[136] = Vec3(0.008000,-0.536000,0.200000);
positions[137] = Vec3(0.088000,-0.061000,0.927000);
positions[138] = Vec3(0.504000,-0.294000,0.910000);
positions[139] = Vec3(-0.860000,0.796000,-0.624000);
positions[140] = Vec3(0.040000,0.544000,-0.748000);
positions[141] = Vec3(0.189000,0.520000,-0.140000);
positions[142] = Vec3(-0.493000,-0.912000,-0.202000);
positions[143] = Vec3(0.815000,0.572000,0.325000);
positions[144] = Vec3(-0.205000,0.604000,-0.656000);
positions[145] = Vec3(0.252000,-0.298000,-0.118000);
positions[146] = Vec3(0.671000,0.464000,-0.593000);
positions[147] = Vec3(0.930000,-0.184000,-0.397000);
positions[148] = Vec3(0.473000,0.500000,0.191000);
positions[149] = Vec3(0.159000,-0.725000,-0.396000);
positions[150] = Vec3(-0.515000,-0.803000,-0.628000);
positions[151] = Vec3(-0.560000,0.855000,0.309000);
positions[152] = Vec3(-0.103000,-0.115000,-0.708000);
positions[153] = Vec3(-0.610000,-0.131000,-0.734000);
positions[154] = Vec3(0.083000,-0.604000,-0.840000);
positions[155] = Vec3(0.688000,-0.200000,-0.146000);
positions[156] = Vec3(0.903000,0.086000,0.133000);
positions[157] = Vec3(-0.136000,0.135000,0.523000);
positions[158] = Vec3(-0.474000,-0.289000,0.477000);
positions[159] = Vec3(0.130000,-0.068000,-0.011000);
positions[160] = Vec3(-0.582000,0.927000,0.672000);
positions[161] = Vec3(0.830000,-0.589000,-0.440000);
positions[162] = Vec3(0.672000,-0.246000,0.154000);
positions[163] = Vec3(-0.212000,-0.142000,-0.468000);
positions[164] = Vec3(-0.021000,0.175000,-0.899000);
positions[165] = Vec3(0.263000,0.326000,0.720000);
positions[166] = Vec3(-0.668000,-0.250000,0.031000);
positions[167] = Vec3(0.822000,-0.860000,-0.490000);
positions[168] = Vec3(0.916000,0.910000,0.291000);
positions[169] = Vec3(-0.358000,-0.255000,0.044000);
positions[170] = Vec3(0.372000,-0.574000,-0.372000);
positions[171] = Vec3(-0.248000,-0.570000,-0.573000);
positions[172] = Vec3(-0.823000,-0.764000,0.696000);
positions[173] = Vec3(-0.848000,0.236000,-0.891000);
positions[174] = Vec3(0.590000,-0.375000,0.491000);
positions[175] = Vec3(-0.153000,0.385000,-0.481000);
positions[176] = Vec3(0.255000,-0.514000,0.290000);
positions[177] = Vec3(0.105000,-0.849000,-0.136000);
positions[178] = Vec3(0.672000,0.203000,-0.373000);
positions[179] = Vec3(0.075000,0.345000,0.033000);
positions[180] = Vec3(-0.422000,0.856000,-0.464000);
positions[181] = Vec3(0.072000,0.166000,0.318000);
positions[182] = Vec3(-0.679000,-0.527000,0.119000);
positions[183] = Vec3(0.613000,0.842000,-0.431000);
positions[184] = Vec3(-0.369000,-0.095000,-0.903000);
positions[185] = Vec3(0.716000,0.565000,-0.154000);
positions[186] = Vec3(-0.412000,-0.642000,-0.229000);
positions[187] = Vec3(0.390000,-0.121000,-0.302000);
positions[188] = Vec3(-0.188000,0.883000,-0.608000);
positions[189] = Vec3(-0.637000,0.325000,0.449000);
positions[190] = Vec3(0.594000,0.745000,0.652000);
positions[191] = Vec3(-0.085000,0.342000,-0.220000);
positions[192] = Vec3(-0.132000,-0.928000,-0.345000);
positions[193] = Vec3(0.859000,-0.488000,0.016000);
positions[194] = Vec3(0.661000,-0.072000,-0.909000);
positions[195] = Vec3(-0.454000,-0.011000,-0.142000);
positions[196] = Vec3(0.859000,-0.906000,0.861000);
positions[197] = Vec3(-0.779000,-0.878000,0.087000);
positions[198] = Vec3(-0.001000,-0.293000,0.851000);
positions[199] = Vec3(0.221000,-0.548000,-0.018000);
positions[200] = Vec3(0.079000,-0.622000,0.653000);
positions[201] = Vec3(0.672000,-0.471000,-0.238000);
positions[202] = Vec3(-0.038000,0.192000,-0.635000);
positions[203] = Vec3(0.428000,0.424000,0.520000);
positions[204] = Vec3(-0.157000,-0.375000,-0.758000);
positions[205] = Vec3(0.317000,0.547000,-0.582000);
positions[206] = Vec3(0.812000,-0.276000,0.687000);
positions[207] = Vec3(-0.438000,0.214000,-0.750000);
positions[208] = Vec3(-0.861000,0.034000,-0.708000);
positions[209] = Vec3(0.770000,-0.532000,0.301000);
positions[210] = Vec3(0.618000,-0.295000,-0.578000);
positions[211] = Vec3(-0.510000,0.052000,0.168000);
positions[212] = Vec3(-0.562000,0.453000,0.691000);
positions[213] = Vec3(-0.269000,0.221000,0.882000);
positions[214] = Vec3(0.039000,-0.785000,0.300000);
positions[215] = Vec3(0.875000,-0.216000,0.337000);
ASSERT_EQUAL_VEC(Vec3( 8.20831e+02, -2.39265e+03, -6.09672e+03), forces1[0], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.66890e+03, 6.12680e+02, 3.26391e+03), forces1[1], TOL);
ASSERT_EQUAL_VEC(Vec3( 8.42139e+02, -5.74020e+03, 9.84681e+02), forces1[2], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.69160e+03, 1.01636e+03, 6.75540e+03), forces1[3], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.04482e+03, 3.89997e+03, 2.79360e+03), forces1[4], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.26801e+03, -3.29937e+03, 6.58525e+03), forces1[5], TOL);
ASSERT_EQUAL_VEC(Vec3( 3.93732e+02, 2.25275e+03, 1.91836e+03), forces1[6], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.61305e+03, -3.55422e+03, -1.60219e+03), forces1[7], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.67011e+03, -5.46823e+03, 2.25854e+03), forces1[8], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.62573e+02, -6.23724e+02, -5.13133e+02), forces1[9], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.46776e+03, -3.41003e+03, 2.99350e+03), forces1[10], TOL);
ASSERT_EQUAL_VEC(Vec3( 7.19389e+02, 4.84737e+03, -3.44817e+03), forces1[11], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.37720e+03, 2.57617e+03, 4.45894e+03), forces1[12], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.38896e+03, 6.92339e+02, -2.28233e+03), forces1[13], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.34724e+03, 1.55928e+03, -8.39021e+02), forces1[14], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.87328e+03, -6.05077e+03, -5.52405e+01), forces1[15], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.59340e+02, 2.11984e+03, 1.90096e+03), forces1[16], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.68172e+03, 3.02879e+03, 3.29144e+03), forces1[17], TOL);
ASSERT_EQUAL_VEC(Vec3(-7.66130e+02, -8.30274e+02, -1.03855e+03), forces1[18], TOL);
ASSERT_EQUAL_VEC(Vec3(-9.02251e+02, 7.56454e+02, 5.31355e+03), forces1[19], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.34128e+03, -2.94917e+03, 3.50106e+03), forces1[20], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.77098e+02, 5.38413e+02, -4.41910e+03), forces1[21], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.89246e+03, -5.97683e+02, -3.30100e+03), forces1[22], TOL);
ASSERT_EQUAL_VEC(Vec3( 9.40266e+02, -4.34341e+03, 4.12621e+02), forces1[23], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.27412e+02, 4.81821e+03, -1.68854e+03), forces1[24], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.63257e+03, -1.35266e+03, -6.23284e+03), forces1[25], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.65374e+03, -1.71701e+03, -4.12940e+02), forces1[26], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.63783e+03, 7.89849e+02, -4.70729e+03), forces1[27], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.95083e+02, 2.48159e+03, 3.71529e+02), forces1[28], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.11286e+03, 4.76531e+03, -1.48722e+03), forces1[29], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.90111e+03, -4.15823e+03, 4.91827e+02), forces1[30], TOL);
ASSERT_EQUAL_VEC(Vec3( 5.31981e+02, -4.61129e+03, -2.26205e+03), forces1[31], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.03402e+03, 2.12277e+03, 4.91126e+03), forces1[32], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.44507e+03, -3.49316e+03, -5.49954e+03), forces1[33], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.39561e+03, 2.02390e+03, 2.72189e+03), forces1[34], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.97484e+02, 1.13306e+03, -1.93722e+03), forces1[35], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.45524e+03, 3.98770e+03, 7.60267e+03), forces1[36], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.18423e+03, 4.30156e+03, 3.13845e+03), forces1[37], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.04712e+03, -4.29254e+03, -1.92390e+03), forces1[38], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.99421e+03, -3.21174e+03, 1.17411e+03), forces1[39], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.00560e+03, 1.25694e+03, -4.55931e+03), forces1[40], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.89329e+03, -1.51189e+01, -1.18920e+03), forces1[41], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.03033e+02, -6.94191e+02, -4.34602e+02), forces1[42], TOL);
ASSERT_EQUAL_VEC(Vec3(-6.55058e+02, 1.86989e+03, -3.20756e+03), forces1[43], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.99529e+03, -8.56961e+02, 2.77210e+02), forces1[44], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.55907e+03, 9.98209e+02, -3.44569e+03), forces1[45], TOL);
ASSERT_EQUAL_VEC(Vec3( 3.15199e+03, -1.77071e+03, 5.33110e+03), forces1[46], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.26890e+03, -2.94837e+03, 8.33300e+02), forces1[47], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.51787e+03, -7.54147e+01, -2.54823e+03), forces1[48], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.76665e+02, 5.52051e+01, -4.82452e+03), forces1[49], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.83919e+03, -3.19345e+03, -1.38614e+02), forces1[50], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.64733e+03, -4.85988e+02, 1.58748e+03), forces1[51], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.83833e+03, 4.19049e+02, -3.08232e+03), forces1[52], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.98287e+03, -3.79608e+02, 2.27447e+03), forces1[53], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.78494e+03, 2.53005e+03, 5.27322e+03), forces1[54], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.31841e+03, 2.35268e+03, 1.34939e+03), forces1[55], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.74641e+03, 2.10630e+03, -6.12271e+03), forces1[56], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.88262e+03, -6.78448e+02, -5.86008e+02), forces1[57], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.78310e+03, 3.67569e+03, -5.48426e+03), forces1[58], TOL);
ASSERT_EQUAL_VEC(Vec3(-5.77653e+02, -7.87374e+02, -6.40935e+03), forces1[59], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.37566e+03, 5.52391e+01, -1.70338e+03), forces1[60], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.99837e+03, 1.04005e+03, -3.77353e+03), forces1[61], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.43485e+03, -3.94193e+03, -3.55960e+03), forces1[62], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.55209e+03, 4.32310e+03, -3.81592e+03), forces1[63], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.46103e+03, -2.98437e+01, 1.00764e+03), forces1[64], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.53255e+03, 5.72328e+02, 3.29519e+03), forces1[65], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.82686e+03, -3.71524e+03, 1.39122e+03), forces1[66], TOL);
ASSERT_EQUAL_VEC(Vec3( 9.04332e+02, -1.44800e+03, 1.57545e+03), forces1[67], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.65856e+03, -2.91481e+03, -6.58907e+03), forces1[68], TOL);
ASSERT_EQUAL_VEC(Vec3( 7.91644e+01, -2.15682e+03, -3.07897e+03), forces1[69], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.95388e+03, 2.61355e+03, 4.74935e+03), forces1[70], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.29404e+03, 3.50134e+03, -4.86504e+02), forces1[71], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.14843e+03, 1.61997e+03, 1.29062e+02), forces1[72], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.61669e+02, 2.47980e+03, 4.78886e+03), forces1[73], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.99354e+03, -1.69720e+02, 2.49340e+03), forces1[74], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.08445e+03, -1.01329e+03, -5.61532e+03), forces1[75], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.35395e+03, -1.97243e+01, -5.03609e+03), forces1[76], TOL);
ASSERT_EQUAL_VEC(Vec3( 3.47258e+03, 1.21225e+03, -2.47694e+03), forces1[77], TOL);
ASSERT_EQUAL_VEC(Vec3( 8.53025e+02, 1.39626e+03, -2.30775e+03), forces1[78], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.74282e+02, 2.21205e+01, 3.40311e+03), forces1[79], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.36776e+02, -4.10780e+01, 2.79555e+03), forces1[80], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.69213e+03, 6.43038e+02, -3.89146e+02), forces1[81], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.01650e+03, -1.17165e+03, -6.24864e+02), forces1[82], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.65143e+03, -1.05607e+03, -2.69230e+03), forces1[83], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.54149e+03, 1.23378e+03, 4.47771e+03), forces1[84], TOL);
ASSERT_EQUAL_VEC(Vec3(-8.11638e+01, 2.00025e+03, -1.71647e+02), forces1[85], TOL);
ASSERT_EQUAL_VEC(Vec3( 9.38381e+02, 1.50858e+03, 3.47900e+03), forces1[86], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.28386e+03, -2.22518e+03, 4.11658e+03), forces1[87], TOL);
ASSERT_EQUAL_VEC(Vec3( 6.58540e+02, 6.72073e+00, 5.59222e+03), forces1[88], TOL);
ASSERT_EQUAL_VEC(Vec3(-7.37375e+02, 1.58093e+03, -4.69233e+03), forces1[89], TOL);
ASSERT_EQUAL_VEC(Vec3( 3.31447e+02, 2.70054e+03, -3.14237e+03), forces1[90], TOL);
ASSERT_EQUAL_VEC(Vec3( 3.57743e+03, 2.12445e+03, 1.15150e+03), forces1[91], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.00033e+03, 8.53614e+02, 1.35134e+03), forces1[92], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.46523e+03, 4.12002e+03, -1.15956e+03), forces1[93], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.47726e+03, -2.12844e+03, -2.30722e+03), forces1[94], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.51392e+03, 2.47607e+02, 1.92133e+03), forces1[95], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.91466e+03, 8.10118e+02, -1.02985e+03), forces1[96], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.37128e+03, -2.45940e+03, 3.70728e+02), forces1[97], TOL);
ASSERT_EQUAL_VEC(Vec3(-5.65127e+01, -3.42531e+01, 9.60430e+03), forces1[98], TOL);
ASSERT_EQUAL_VEC(Vec3( 6.15217e+02, -1.68773e+03, 1.42295e+03), forces1[99], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.09417e+02, 2.10263e+03, 1.04784e+03), forces1[100], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.93317e+03, 2.06687e+03, -1.64072e+03), forces1[101], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.36623e+03, -2.43368e+03, 1.57335e+03), forces1[102], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.66360e+03, -3.56116e+03, 4.66828e+02), forces1[103], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.46239e+03, -2.54761e+03, 7.26200e+02), forces1[104], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.48952e+03, 2.17361e+03, -2.35016e+03), forces1[105], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.04960e+03, 2.83274e+03, 4.20052e+03), forces1[106], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.08299e+03, 1.18791e+03, -4.67708e+03), forces1[107], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.42314e+03, -2.95658e+03, -1.31832e+03), forces1[108], TOL);
ASSERT_EQUAL_VEC(Vec3( 3.76388e+03, 4.33363e+03, 3.15347e+03), forces1[109], TOL);
ASSERT_EQUAL_VEC(Vec3(-8.49993e+02, 1.28674e+03, -5.50558e+03), forces1[110], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.77761e+02, -1.41841e+03, 2.92791e+02), forces1[111], TOL);
ASSERT_EQUAL_VEC(Vec3(-7.08175e+02, -5.37887e+03, 3.78340e+03), forces1[112], TOL);
ASSERT_EQUAL_VEC(Vec3(-5.91441e+02, 2.10098e+03, -4.82150e+00), forces1[113], TOL);
ASSERT_EQUAL_VEC(Vec3(-6.92132e+02, 4.62800e+03, 4.45194e+03), forces1[114], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.31335e+03, 3.51626e+03, -2.57372e+03), forces1[115], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.55816e+03, -1.01211e+03, -5.58678e+03), forces1[116], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.24174e+03, 1.25094e+03, -6.95682e+02), forces1[117], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.93053e+03, -1.12989e+03, -1.22677e+03), forces1[118], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.86217e+02, 4.00249e+02, -2.23733e+03), forces1[119], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.26876e+03, -3.20779e+03, -4.08934e+03), forces1[120], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.18596e+03, -7.37298e+02, 7.15904e+03), forces1[121], TOL);
ASSERT_EQUAL_VEC(Vec3( 3.95098e+03, -3.37054e+03, -2.70316e+03), forces1[122], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.35471e+03, 2.40580e+03, 1.84528e+03), forces1[123], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.61347e+03, -3.36020e+03, -3.99614e+03), forces1[124], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.72566e+03, -1.88479e+03, -1.69932e+03), forces1[125], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.62207e+03, 1.43319e+03, 7.47053e+02), forces1[126], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.80955e+03, -4.54150e+02, 7.21597e+03), forces1[127], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.56885e+01, 1.48364e+03, -1.62205e+03), forces1[128], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.47403e+02, -2.07399e+03, 1.49903e+03), forces1[129], TOL);
ASSERT_EQUAL_VEC(Vec3( 9.86979e+01, 2.60134e+03, 2.12407e+02), forces1[130], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.27811e+03, 3.71499e+03, -4.95036e+03), forces1[131], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.80783e+03, -1.54371e+03, -4.77208e+03), forces1[132], TOL);
ASSERT_EQUAL_VEC(Vec3(-9.61648e+01, -2.59415e+03, 1.96258e+03), forces1[133], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.89381e+02, 4.81189e+03, 9.17885e+02), forces1[134], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.93800e+02, 1.38919e+03, -2.46053e+03), forces1[135], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.62694e+03, 2.35155e+03, 3.91729e+03), forces1[136], TOL);
ASSERT_EQUAL_VEC(Vec3( 5.31966e+03, -8.44013e+02, -3.54397e+03), forces1[137], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.35304e+02, -3.23988e+03, -2.69527e+03), forces1[138], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.44149e+03, -2.15586e+03, 3.15946e+03), forces1[139], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.23043e+02, 2.28842e+03, -3.86313e+03), forces1[140], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.70608e+02, 7.17360e+02, 1.22594e+03), forces1[141], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.21938e+03, 8.98142e+02, 1.48998e+03), forces1[142], TOL);
ASSERT_EQUAL_VEC(Vec3(-6.45287e+02, -9.68214e+02, 1.22919e+03), forces1[143], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.76443e+03, 2.83867e+02, -4.18068e+03), forces1[144], TOL);
ASSERT_EQUAL_VEC(Vec3(-6.06308e+02, 1.35580e+03, 3.57015e+03), forces1[145], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.13396e+03, -2.12738e+03, 4.47039e+00), forces1[146], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.45064e+03, 7.41311e+02, -4.13913e+03), forces1[147], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.08512e+02, 1.56784e+03, 4.52361e+03), forces1[148], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.01824e+02, -2.64844e+03, -2.53109e+03), forces1[149], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.89630e+03, 4.60390e+03, -2.34064e+03), forces1[150], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.39723e+03, 1.74535e+03, 3.47741e+03), forces1[151], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.18808e+03, -1.14056e+03, 5.07060e+02), forces1[152], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.47101e+03, -2.10613e+03, -1.86346e+03), forces1[153], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.73424e+03, -4.59311e+03, -6.97340e+03), forces1[154], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.25831e+02, 2.12782e+03, 6.03112e+02), forces1[155], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.55032e+03, 3.85747e+03, 5.93351e+03), forces1[156], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.28816e+02, -3.00207e+03, 2.10094e+02), forces1[157], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.23671e+03, -1.00832e+03, 2.76637e+03), forces1[158], TOL);
ASSERT_EQUAL_VEC(Vec3(-6.46593e+02, 1.33324e+03, 4.52795e+03), forces1[159], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.37177e+03, 2.53533e+02, -7.78677e+02), forces1[160], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.47620e+03, -1.99006e+02, -3.16482e+03), forces1[161], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.00679e+03, 1.70898e+03, 4.88103e+03), forces1[162], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.52761e+03, -5.32723e+02, 4.40420e+03), forces1[163], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.92570e+03, 3.00581e+03, -4.10543e+03), forces1[164], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.62169e+02, 1.69864e+03, -8.35330e+02), forces1[165], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.40942e+03, -5.87402e+01, 2.05229e+03), forces1[166], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.77206e+03, 2.68955e+02, 9.28640e+02), forces1[167], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.56436e+03, -8.49851e+02, 2.03495e+03), forces1[168], TOL);
ASSERT_EQUAL_VEC(Vec3( 3.17554e+03, -2.75544e+03, 3.21077e+03), forces1[169], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.85777e+03, -1.86570e+03, 1.68128e+02), forces1[170], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.03897e+03, 1.61130e+02, -1.77268e+03), forces1[171], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.72766e+02, 1.73951e+03, -1.87733e+03), forces1[172], TOL);
ASSERT_EQUAL_VEC(Vec3(-6.55527e+02, 2.24951e+03, -4.12361e+03), forces1[173], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.02507e+03, -1.03114e+03, 3.14422e+03), forces1[174], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.43834e+03, 3.75451e+02, 2.19531e+03), forces1[175], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.78662e+03, 3.61832e+02, 6.06302e+03), forces1[176], TOL);
ASSERT_EQUAL_VEC(Vec3( 3.01198e+02, -1.53666e+03, 2.85060e+03), forces1[177], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.07398e+03, 6.61994e+02, -8.61972e+02), forces1[178], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.54259e+03, 9.35531e+02, 2.58246e+03), forces1[179], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.22576e+03, -3.02376e+03, -7.71072e+02), forces1[180], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.76413e+03, -2.59918e+03, 9.81425e+02), forces1[181], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.99912e+03, -2.08518e+03, 2.96759e+03), forces1[182], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.05925e+03, -7.79905e+02, -6.83796e+01), forces1[183], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.17328e+03, -3.35387e+03, -5.62594e+03), forces1[184], TOL);
ASSERT_EQUAL_VEC(Vec3( 4.34398e+03, 1.27413e+03, 1.09523e+03), forces1[185], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.38740e+03, 3.55672e+03, 1.53500e+03), forces1[186], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.14439e+03, 3.47927e+03, -2.29435e+03), forces1[187], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.52579e+03, 1.41487e+03, -6.37862e+03), forces1[188], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.17319e+03, -3.07753e+03, 2.16801e+03), forces1[189], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.61247e+03, 1.39162e+03, -3.30074e+03), forces1[190], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.53680e+03, -5.26042e+02, 1.85855e+03), forces1[191], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.88755e+03, 3.66622e+02, 1.40663e+02), forces1[192], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.97756e+02, -5.48803e+03, 9.38250e+02), forces1[193], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.26511e+03, 1.95703e+03, -3.01994e+03), forces1[194], TOL);
ASSERT_EQUAL_VEC(Vec3( 8.24840e+02, 3.33922e+02, -2.23925e+03), forces1[195], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.16139e+03, 3.33085e+03, -2.67863e+03), forces1[196], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.94535e+03, 4.31827e+02, -1.43010e+03), forces1[197], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.15722e+02, -3.67489e+03, -6.60062e+03), forces1[198], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.08549e+03, -1.10080e+03, 3.87260e+03), forces1[199], TOL);
ASSERT_EQUAL_VEC(Vec3(-7.34195e+02, -2.62191e+03, -6.68007e+02), forces1[200], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.43085e+02, -4.49257e+03, 1.69265e+03), forces1[201], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.61784e+03, 1.36449e+02, 1.89887e+03), forces1[202], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.77060e+03, 1.21537e+03, -1.61977e+03), forces1[203], TOL);
ASSERT_EQUAL_VEC(Vec3(-3.11482e+03, -3.10409e+03, -3.86688e+03), forces1[204], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.21598e+02, 5.26726e+02, 1.90196e+03), forces1[205], TOL);
ASSERT_EQUAL_VEC(Vec3( 1.86923e+03, 3.90716e+02, 1.66422e+03), forces1[206], TOL);
ASSERT_EQUAL_VEC(Vec3(-4.87140e+03, 4.24474e+03, -6.93148e+02), forces1[207], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.07292e+03, -6.60370e+01, -1.35550e+03), forces1[208], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.62135e+03, -4.27394e+03, 3.63843e+03), forces1[209], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.91469e+03, 1.89371e+02, -4.65934e+03), forces1[210], TOL);
ASSERT_EQUAL_VEC(Vec3(-7.40876e+02, -1.17545e+03, 5.03240e+03), forces1[211], TOL);
ASSERT_EQUAL_VEC(Vec3(-5.36901e+02, 1.95912e+03, 2.90388e+03), forces1[212], TOL);
ASSERT_EQUAL_VEC(Vec3(-2.55522e+03, 4.34339e+03, -5.12296e+03), forces1[213], TOL);
ASSERT_EQUAL_VEC(Vec3(-1.08022e+03, -4.69890e+03, 6.58358e+03), forces1[214], TOL);
ASSERT_EQUAL_VEC(Vec3( 2.42115e+03, 2.85514e+03, 5.90927e+03), forces1[215], TOL);
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment